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ABSTRACT 

The  objective  of  this  thesis  is  to  test  the  Semianalytical  Satellite 
Theory  (SST)  as  implemented  in  The  Charles  Stark  Draper  Laboratory  version 
of  the  Goddard  Trajectory  Determination  System  (GTDS)  against  long  arcs  of 
real  data  for  highly  eccentric  orbits.  The  real  data  is  in  the  form  of 
North  American  Defense  Command  (NORAD)  element  sets  and  actual  observa¬ 
tions.  Data  are  pre-processed  and  converted  to  a  set  of  observations  in  a 
GTDS -compatible  format.  These  observations  undergo  a  differential  correc¬ 
tion  (DC)  process  to  generate  the  initial  conditions  for  an  SST  ephemeris 
prediction.  This  prediction  is  then  compared  with  real  data  to  evaluate 
the  performance  of  the  SST. 


Data  available  for  this  study  included  element  sets  and  metric  observa¬ 
tions  for  the  following  objects: 


Object 


Data  Type 


Period 


NSSC  9829 
Molniya 

( 1 977-1 0A) 

2-17 

Mean  elements 

February 

1977  -  July 

1986 

NSSC  14095 

Exosat 

( 1983-51A) 

Mean  elements 

Observations 

June  1983 
June  1983 

-  December 

-  December 

1985 

1985 

NSSC  13964 

( 1  983-25A ) 

Mean  elements 

May  1983 

-  September 

1985 

Molniya 

1-57 

Observations 

May  1983 

-  December 

1985 

The  results  validated  the  performance  of  the  Semianalytical  Satellite 
Theory  for  high  eccentricity  orbits.  When  elements  sets  were  used  as 
inputs  to  the  DC,  comparisons  between  the  NORAD  singly-averaged  elements 
and  the  SST  predictions  showed  good  agreement.  When  the  metric 
observations  were  used  as  inputs  to  the  DC,  agreement  between  the  SST 
predictions  and  the  observation  data  was  very  good. 
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CHAPTER  1 


INTRODUCTION 

The  objective  of  this  thesis  is  to  test  the  Semianalytical  Satellite 
Theory  (SST)  as  implemented  in  the  Charles  Stark  Draper  Laboratory  version 
of  the  Goddard  Trajectory  Determination  System  (GTDS)  against  long  arcs  of 
real  data  for  highly  eccentric  orbits.  The  real  data  is  in  the  form  of 
North  American  Defense  Command  (NORAD)  element  sets  and  metric  observa¬ 
tions.  The  orbits  of  primary  interest  will  be  those  of  two  Soviet  Molniya 
spacecraft  (12-hrf  high  e)  and  of  the  European  Space  Agency  satellite 
Exosat.  Data  are  pre-processed  and  converted  to  a  set  of  observations  in 
a  GTDS -compatible  format.  These  observations  undergo  a  differential 
correction  (DC)  process  to  generate  the  initial  conditions  for  an  SST 
ephemeris  prediction.  Processing  of  the  real  data  proceeds  along  slightly 
different  paths  depending  on  whether  element  sets  or  observations  are 
used  to  generate  the  initial  conditions.  This  prediction  is  then  compared 
with  real  data  to  evaluate  the  performance  of  the  SST. 

Chapter  1  reviews  the  use  of  elliptical  orbits  by  the  scientific  and 
military  communities.  The  need  for  long-term  orbit  predictions  is 
described  and  the  application  of  artificial  satellite  theory  to  highly 
elliptical  orbits  is  reviewed. 


test irq  new  instruments  or  engineering  techniques,  and  a  few  have  been  for 
the  mapping  ;t  Earth  resources  [1]. 

The  primary  uses  of  elliptical  orbits  have  been  for  scientific  missions 
-oniutel  py  all  the  space-going  countries  of  the  world,  military  missions 
conducted  p:  mar;  iy  by  the  United  States  and  the  Soviet  Union,  and  for 
transfer  orbi's  r  geosynchronous  and  other  high  altitude  orbits.  The 
ma  •  r  advantages  of  highly  elliptical  orbits  include  relatively  long 

“s  *  i . -?  .1  relatively  small  area  of  space  and  an  indep  idence  from 

a--  sp'.er.  Ir-au  ef !«  -fs  during  most  of  the  orbit  lifetime.  These  advan- 
*  -.av*  he,"  ex;  :•  ;  ted.  by  both  scientific  and  military  investigators. 

lcier.tif:  -  Mi  s  s  i  ins 

T-e  jrea  *  ma '  •*.  r  i  ♦  y  of  scientific  satellites  operated  in  elliptical 
n;  v;:  ns  h.  ave  bee;,  devote  1  to  the  studies  of  the  Earth  and  the  near-Earth 
er.y.  r  .-rimer  •  .  The  primary  users  have  been  the  United  States,  the  Soviet 
’:M  **-•*.,  aril  the  European  Space  Agency  (ESA). 

n»-  it'  th"  early  scientific  programs  was  the  Interplanetary  Monitor 
E  latf  rrr  (  IMl  program  which  consisted  of  ten  United  States  spacecraft 

launched  fr  >m  1  'if- 1  to  'W7 .  Its  mission  was  to  continuously  measure  the 
rail i’ i  it  environment  in  the  immediate  vicinity  of  Earth  and  in  inter¬ 

planetary  spa  'e  during  a  complete  solar  cycle  in  preparation  for  the 
Ap  11  ir  iram.  The  satellite  orbits  had  a  typical  perigee  height  of 

apt  r  x  ;  y  2  k  t  and  an  apogee  height  of  approximately  200,000  km  [2]. 


A  more  current  satellite,  the  Combined  Release  and  Radiation  Effects 
Satellite  (CRRES)  has  an  unusual  mission  in  that  it  performs  related 
scientific  experiments  in  both  low  earth  orbit  and  a  geosynchronous 
transfer  orbit  [3].  One  mission  objective  of  the  CRRES  is  to  determine 
how  the  Earth's  magnetic  field  affects  in  situ  microelectronic  components. 
Another  objective  is  to  gather  the  data  for  an  improved  model  of  the 
Earth's  radiation  environment.  The  orbital  parameters  of  the  three  year 
mission  include  a  geosynchronous  altitude  apogee,  a  low  inclination,  and  a 
perigee  altitude  of  approximately  400  km. 

The  Soviet  Union  continues  its  active  use  of  space,  and  several  recent 
spacecraft  have  utilized  the  elliptical  orbit.  The  Astron  satellite 
launched  in  March  1983  from  Tyuratam  had  a  mission  to  investigate  distant 
X-ray  and  ultraviolet  emissions.  Its  orbit  measured  2000  km  by  200,000  km 
and  had  an  orbital  period  of  four  days.  It  thus  remained  at  great 
distances  from  the  Earth  for  three  and  a  half  days  of  its  orbit  and  was 
highly  visible  to  Soviet  ground-based  tracking  stations.  The  Prognoz 
series  began  in  1972  and  has  included  a  wide  range  of  scientific  missions. 
Prognoz  9  was  launched  in  July  1983  from  Tyuratam  into  an  orbit  measuring 
400  km  by  720,000  km  to  investigate  remnant  radiation  from  the  "big  bang". 
Prognoz  10  was  launched  in  April  1985  into  an  orbit  measuring  400  km  by 
200,000  km  to  study  the  shock  wave  boundaries  created  by  the  interaction 
of  the  solar  wind  plasma  and  the  Earth's  magnetosphere.  Other  Soviet 
scientific  satellite  programs  planned  for  later  in  the  decade  [4]  include 
placing  the  first  of  two  pair  of  Prognoz-type  and  Magion-type  satellites 
into  an  orbit  measuring  500  km  by  250,000  km  and  the  second  pair  into  an 
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orbit  with  an  apogee  of  only  20,000  km.  The  first  set  will  investigate 
the  magnetospheric  tail  as  the  other  pair  observes  the  auroral  field 
lines . 

The  ESA  began  utilizing  the  elliptical  orbit  in  August  1975  with  its  Cos-B 
spacecraft  designed  as  a  gamma-ray  observatory.  It  initially  operated  in 
a  0.88  eccentricity  orbit  of  340  km  by  99,870  km  and  decayed  in  January 
1986.  The  ISEE-2,  launched  in  October  1977,  was  designed  for  magneto- 
spheric  exploration  from  a  0.91  eccentricity  orbit.  It  is  predicted  to 
decay  in  September  1987  [5].  Another  ESA  satellite,  whose  orbit  will  be 
studied  in  detail  later  in  this  paper,  is  the  Exosat,  launched  in  May 
1983. 

1.1.2  Military  Missions 

The  traditional  military  missions  in  space  are:  (1)  communications;  (2) 
reconnaissance  and  surveillance;  (3)  navigation;  (4)  meteorology;  and  (5) 
geodesy.  Spacecraft  with  these  missions  primarily  operate  in  one  of  five 
orbit  types:  (1)  low  Earth  orbits,  (LEO)  below  5000  km;  (2)  geosynchronous 
orbits  (GEO),  at  35,700  km  altitude;  (3)  Molniya  orbits,  about  500  km  by 
40,000  km  with  a  12  hour  period;  (4)  semi -synchronous  orbits,  at  20,000  km 
with  a  12  hour  period;  and  (5)  super-synchronous  orbits  between  GEO  and 
the  Moon.  Tables  1  thru  4  of  [6]  describe  the  current  U.S.  and  Soviet 
military  space  systems.  The  U.S.  utilizes  elliptical  orbits  for  its 
Satellite  Data  System  (SDS),  and  the  Soviets  utilize  elliptical  orbits  for 
their  Molniya  communications  satellites  as  well  as  some  early  warning 
missions . 


There  is  little  open  literature  on  the  classified  SDS  program.  Reference 
(7]  describes  it  as  having  a  three-fold  mission:  (1)  maintain  polar 
communication  with  Strategic  Air  Command  aircraft  in  those  areas  uncovered 
by  communication  satellites  in  geosynchronous  orbit;  (2)  provide  a  command 
and  control  link  for  USAF  satellites;  and  (3)  act  as  a  relay  for  KH-1 1 
photographic-reconnaissance  satellites.  The  SDS  spacecraft  are  launched 
atop  the  Titan  3B  /  Agena  booster  and  inserted  into  an  orbit  approximately 
320  km  by  38,600  km,  with  apogee  timed  to  occur  over  the  polar  regions. 
An  inclination  of  about  63°  provides  greater  satellite  visibility  at  high 
latitudes.  There  were  an  estimated  eight  operational  spacecraft  launched 
between  1975  and  1981  [7], 

The  Soviet  Union  inaugurated  its  eminently  successful  domestic  communica¬ 
tion  satellite  system  in  April  of  1965  with  the  launch  of  Molniya  1-1. 
Since  then  1 1 0  of  the  1600  kg  spacecraft  have  carried  the  Molniya  name. 

Due  to  the  enormous  expanse  of  the  Soviet  Union  -  over  18C°  in  longitude 
and  over  40°  in  latitude  to  the  Arctic  Circle  -  the  Soviets  pioneered  the 
use  of  highly  elliptical  orbits  inclined  63°  to  66”  to  the  equat  >r  to 
provide  reliable,  nation-wide  telephone,  telegraph,  and  television 
communications.  With  perigee  heights  of  400  to  600  Km  situated  over  the 
Southern  Hemisphere,  Molniya  satellites  linger  over  the  Soviet  Union  at 
altitudes  of  39,000  to  40,000  km  for  approximately  eight  hours  every 
revolution.  By  carefully  positioning  satellites  in  sequence  as  few  as 
three  satellites  can  provide  around-the-clock  coverage,  although  four  or 
more  are  usually  employed. 


Two  constellations  of  Molniya  satellites  serve  the  Soviet  communications 
requirements  [4].  The  older  Molniya  1  satellites  are  flown  in  eiqht 
orbital  planes  spaced  45°  apart.  Each  satellite  is  positioned  in  its 
plane  to  ensure  that  every  satellite  in  the  network  traces  the  same  path 
over  the  Earth  each  day.  This  technique  minimizes  trackinq  requirements 
of  the  more  than  90  Orbita  qround  stations.  The  second  system  is 
comprised  of  advanced  Molniya  3  satellites,  which  are  thouqht  tr  be 
enhanced  for  maritime  communi cations .  The  Molniya  3  network  presently 
consists  [4]  of  four  satellites  spaced  90°  apart  with,  ascendinq  nodes  near 
bc.i°  and  four  satellites  interspersed  with  ascendinq  nodes  near  'SS0. 
K:  )ur-  1  [4]  is  a  diaqram  of  the  Molniya  3-21,  and  Kiqure  2  [d]  shows  a 
typical  Molniya  qroundt rack  . 
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1 . 1 .3  Transfer  to  Geosynchronous  Orbit 
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than  for  scientific  satellites  because  of  the  increased  influence  of  air 


drag  at  perigee  due  to  their  large  size.  Long-term  predictions  also  can 
be  particularly  difficult  if  the  lunar-solar  perturbations  hold  the 
perigee  height  in  the  atmosphere,  causing  the  satellite  to  pass  through 
several  resonance  conditions. 


One  such  object  in  the  geosynchronous  transfer  orbit  is  NSSC  10723 
( 1 978-1 2C),  third  stage  for  the  International  Ultraviolet  Explorer 
launched  in  January  1978.  It  had  an  initial  semimajor  axis  of  29,648  km 
and  perigee  height  of  177  km,  for  an  eccentricity  of  0.779.  It  had  an 
estimated  orbital  lifetime  of  seven  years  [1]. 


1.2  NEED  FOR  LONG-TERM  PREDICTIONS 


There  are  several  reasons  for  investigating  the  behavior  of  highly 
elliptical  orbits.  The  importance  of  accurately  predicting  the  long-term 
motion  of  a  satellite  orbit  is  probably  most  evident  during  the  mission 
design  phase  of  the  spacecraft.  Besides  the  concerns  that  exist  for  the 
orbit  evolution  of  a  single  spacecraft,  there  are  concerns  for  the  evolu¬ 
tion  over  time  for  a  constellation  of  two  or  more  spacecraft.  The  need  to 
reliably  predict  satellite  orbits  is  evident  in  light  of  the  increasing 
attention  given  to  the  problem  of  orbital  debris  and  to  the  possibility  of 
collisions  in  space.  Long-term  prediction  of  satellite  orbits  has  also 


become  an  issue  in  the  area  of  space  arms  control 


1.2.1  Mission  Design 


The  orbital  lifetime  should  be  well  understood  prior  to  launch  via  the 
mission  analysis  process.  Obviously,  the  satellite  orbital  lifetime  must 
be  longer  than  the  satellite  operational  lifetime.  Long-term  predictions 
of  the  orbit  insure  the  satellite  has  the  proper  orbital  parameters  while 
the  spacecraft  is  performing  its  intended  function.  Long-term  predictions 
are  also  necessary  to  design  the  satellite  stationkeeping  algorithms. 
Large  and  heavy  satellites  do  not  completely  burn  up  during  atmospheric 
re-entr  •  the  exact  time  and  place  of  re-entry  must  be  predicted  in  order 
to  warn  the  inhabitants  of  any  danger  (e.g.,  Cosmos  954  in  1978,  Sky  lab  in 
1979,  and  Cosmos  1402  in  1983). 

The  most  significant  sources  of  perturbations  of  geocentric  satellite 
orbits,  and  especially  for  mission  design  studies,  are  the  central  body 
oblateness,  atmospheric  drag,  and  lunar-solar  gravity.  The  oblateness  of 
the  Earth  causes  a  secular  change  in  orbit  orientation  and  is  important 
directly  because  of  its  immediate  influence  on  the  mission,  and  indirectly 
because  by  changing  the  orbital  relationship  to  the  lunar-solar  constella¬ 
tion,  it  alters  the  lunar-solar  effects.  Atmospheric  drag  results  in  a 
decrease  in  orbital  energy  and  therefore  a  secular  reduction  of  apogee. 
For  orbits  of  higher  ecentricity,  with  e  approximately  equal  to  0.7,  the 
gravitational  perturbations  due  to  the  Sun  and  Moon  are  dominant  in  deter¬ 
mining  the  evolution  of  satellite  orbits.  The  orbital  elements  signifi¬ 
cantly  effected  by  external  body  perturbations  are  the  eccentricity, 
inclination,  argument  of  perigee,  and  the  longitude  of  the  ascending  node. 
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The  perigee  variations  due  to  the  Sun  and  Moon  follow  a  characteristic 


pattern  approximated  by  the  superposition  of  four  sine  waves:  (1)  a  very 
small -amplitude  wave  with  the  same  periodicity  as  the  satellite;  (2)  a 
14-day  wave  of  larger  amplitude;  (3)  a  183-day  wave  of  amplitude  5.8  times 
as  large;  and  (4)  a  several-year  period  wave,  usually  of  still  larger 
amplitude.  The  four  waves  are  usually  called  the  short,  intermediate, 
long,  and  very-long  period  effects  [10].  The  last  two  are  important  for 
determining  orbital  lifetimes  of  high-apogee,  low-perigee  orbits.  The 
intermediate  and  long-period  effects  are  related  to  the  perturbing  body's 
orbital  motion  and  the  very-long  effect  to  the  satellite  orbit  plane's 
changing  orientation. 

1.2.2  Catalog  Maintenance 

The  responsibility  for  space  population  record-keeping  is  entrusted  to  the 
North  American  Aerospace  Defense  Command  (NORAD).  Observations  on  the 
deep  space  satellites  are  provided  by  four  sources  -  the  optical  sensors 
at  the  Baker-Nunn  camera  sites,  electro-optical  sensors  including  the 
three  Ground  Based  Electro-Optical  Deep  Space  Sensor  System  (GEODSS) 
sites,  deep  space  radar  sensors,  and  near-Earth  radars  such  as  NAVSPASUR 
and  PAVE  PAWS  when  the  satellites  pass  with  range  of  these  sensors. 
Observations  from  these  sensors  are  used  in  conjuction  with  high-speed 
computers  for  identification  and  cataloging  of  all  satellites  and  objects 
in  space.  The  resulting  element  sets  are  then  distributed  to  users  to 
allow  prediction  of  the  future  position  and  velocity  of  each  space  object. 
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Deep  space  satellites  present  a  difficult  problem  in  terms  of  catalog 
maintenance  because  many  of  the  satellite  orbits  are  out  of  range  or  at 
the  very  limits  of  present  tracking  capabilities  [11].  Because  of  their 
slow  angular  motion  relative  to  near-Earth  satellites,  deep  space  satel¬ 
lites  are  observed  by  fewer  sensors  and  their  orbits  are  therefore 
calculated  with  less  data. 

Reference  [12]  stated  that  the  catalog  of  objects  in  space  is  relatively 
complete  for  objects  greater  than  10  cm  in  size,  and  he  estimates  that  the 
population  of  objects  1  -  4  cm  in  size  could  be  3  to  10  times  the  known 

population.  An  object  this  size  could  have  a  mass  of  100  grams  and  at 

relative  velocities  of  the  order  of  10  km/sec  could  inflict  significant 
damage  to  a  satellite  or  an  abandoned  rocket  stage.  Reference  [12]  went 

on  to  state  that  it  is  possible  to  detect  objects  1-10  cm  across  using 
space-based  radar  or  lidar  systems,  but  that  the  mass  and  power 

requirements  are  currently  prohibitively  high. 

1.2.3  Satellite  Constellation  Evolution 

The  parameters  of  satellite  constellation  design  include  the  number  of 
satellites  available,  the  properties  of  the  sensor  and  payload  (constrain¬ 
ed  by  the  operating  altitude  and  attitude  of  the  spacecraft),  and  the 
amount  and  type  of  coverage  required.  Methods  of  minimizing  revisit  time 
or  maximizing  coverage  typically  assume  symmetric  constellations  [13]. 
That  is,  satellites  in  a  single  orbit  plane  are  evenly  distributed  in  mean 
anomaly,  and  individual  orbit  planes  at  common  inclinations  are  evenly 
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distributed  in  right  ascension  of  ascending  node.  Real -world  perturba¬ 
tions  can  and  do  cause  changes  in  the  argument  of  perigee,  loncitude  of 
the  ascending  node,  etc.  These  real-world  perturbations  necessitate 
careful  study  of  the  long-term  stability  of  the  satellite  constellation. 

1.2.4  Debris  and  Collision  Assessments 

More  than  3000  payloads  and  12,000  objects  have  been  put  into  orbit  since 
1957.  The  December  1984  "Satellite  Situation  Report",  quoted  in  [12], 
listed  5408  objects  in  orbit,  of  which  5%  were  operational  payloads,  20% 
were  non-operational  payloads,  25%  were  mission-related  objects,  and  50% 
were  satellite  breakups.  Many  mission-related  objects  -  upper  stages, 
apogee  kick  motors,  spent  rocket  casings,  parts  of  separation  devices, 
protective  covers,  etc.  -  remain  in  low  Earth  orbit  or  in  the  geo¬ 

synchronous  transfer  orbit.  The  geos  .ationary  ring,  a  strictly  limited 
territory  only  264,930  km  long,  is  becoming  increasingly  crowded, 

especially  near  the  stable  points  of  107°  West  and  76°  East  longitude. 
Trends  indicating  that  the  number  of  operating  spacecraft  will  double 
every  five  years  [14]  accentuate  the  need  for  accurate  long-term 
prediction  of  objects  in  these  orbits. 

There  are  several  means  of  reducing  the  hazards  of  collision.  Some  are 
relatively  easy  to  implement  with  spacecraft  design  changes  or  operational 
procedures.  These  include  reducing  the  number  of  objects  at  separation, 

de-orbiting  large  satellites  to  re-enter  the  atmosphere,  depleting 

residual  propellants,  and  conducting  satellite  destruction  tests  at  low 
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altitude  to  ensure  rapid  decay  of  debris.  Another  method,  which  is  less 
costly  in  terms  of  the  required  spacecraft  mass  than  deorbiting  to  Earth 
is  to  desynchronize  upward.  It  is  written  policy  of  the  INTELSAT  board  of 
governors  to  desynchronize  unusable  spacecraft  40  km  to  50  km  upward;  it 
is  written  policy  of  NOAA  to  desynchronize  at  least  300  km  upward  [14]. 
The  Soviet  Union  does  not  circularly  desynchronize,  but  has  fired  to 
produce  large  eccentricity  to  reduce  the  chances  of  crossing  on  the 
geostationary  ring. 

1.2.5  Space  Arms  Control  Issues 

The  requirement  for  long-term  prediction  of  orbits  is  relevant  in 
discussions  of  military  space  systems.  The  vulnerability  of  a  satellite 


to  attack 

is  primarily  a  function  of 

the 

orbit 

in 

which 

the  satellite 

operates , 

since  the  orbit  determines 
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cos  t 

for  an  Earth- 

launched  ASAT  intercept  or  the  power  for  a  laser  or  particle  beam  attack. 

There  has  been  general  agreement  that  essential  warning  and  communication 
capabilities  should  be  preserved  at  a  very  high  priority  even  in  the  face 
of  a  nuclear  attack  [6,15,16].  One  possible  bilateral  agreement  to  en¬ 
hance  the  security  of  seme  space  assets  most  relevant  to  nuclear  stability 
in  an  adequately  verifiable  way  could  develop  "rules  of  the  road"  for 
space,  i.e.,  an  international  agreement  governing  the  use  of  space.  This 
agreement  might,  for  example,  make  it  legal  for  any  signatory  nation  to 
take  seme  action  regarding  other  nations'  satellites  coming  closer  than 
certain  stated  distances  to  specified  satellites  of  its  own.  The  stated 


distances  and  specified  satellites  would  have  to  be  resolved  according  to 
the  orbits  of  the  satellites  to  be  protected  and  according  to  tf>- 
requirements  of  civilian  space  users. 

Illustrative  "rules  of  the  road"  include:  (1)  Keep  TOO  Km  and  3°  out-of¬ 
plane  from  foreign  satellites  below  5000  Km;  (2)  Keep  500  Km  from  foreign 
satellites  above  5000  Km  except  those  within  500  Km  of  geosynchronous 
altitude;  (3)  allow  one  pre -announced  close  approach  at  a  time;  (4)  find 
guilty  of  trespass  the  satellite  of  the  nation  which  most  recently 
initiated  a  maneuver  burn  [16]. 

1.3  APPLICATION  OF  ARTIFICIAL  SATELLITE  THEORY  TO  ELLIPTICAL  ORBITS 

The  application  of  analytical,  senianalytical,  and  numerical  orbit  propa¬ 
gators  to  the  problem  of  long-term  prediction  of  eccentric  orbits  is  the 
subject  of  this  section.  When  making  any  comparison,  it  is  important  tc 
understand  the  two  things  being  compared.  Because  one  objective  of  this 
research  was  to  test  the  long-term  predictions  of  the  semianalytical 
theo  y  against  real  data  in  the  form  of  NORAD  mean  elements,  both  the 
semianalytical  and  NORAD  satellite  theories  will  be  reviewed.  Also 
included  is  a  review  of  previous  applications  of  artificial  satellite 
theory  to  the  problem  of  highly  elliptical  orbits. 
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Reference  [2b]  described  the  basic  theory  for  a  computer  program  PROD 
Program  for  Orbit  Development)  which  predicted  the  long-term  development 
•  f  elliptic-  orbits  under  the  influence  of  the  Earth's  gravitational 
p  cenr;  a',  ani  the  attractions  of  the  Sun  and  Moon.  Effects  of  drag  and 
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solar  radiation  pressures  were  not  included.  Reference  [21]  used  this 
program  to  predict  orbital  elements  for  a  Molniya  1  satellite  and  compared 
them  to  U.S.  Air  Force  orbital  5-line  elements.  The  PROD-predicted  peri¬ 
gee  heights  were  accurate  to  about  20  km  throughout  a  21 -month  interval. 

Reference  [21]  also  developed  an  analytical  method  for  1 2-hr  orbits  of  65° 
inclination  and  demonstrated  that  the  lifetime  must  be  between  one  and 

seven  years.  It  was  shown  that  in  the  absence  of  air  drag,  orbital  life¬ 

time  is  controlled  by  the  variation  in  perigee  height  which,  in  turn,  is 
controlled  by  the  variation  in  eccentricity.  This  variation  in  eccen¬ 
tricity  due  to  lunar-solar  perturbations  was  expressed  as  a  function  of 
argument  of  perigee  and  longitude  of  the  ascending  node.  With  initial 

values  of  eccentricity  and  argument  of  perigee  (e=0.74  and  o)=286°)  for  31 
of  the  65  Molniyas  launched  up  to  that  time,  and  with  an  observed  value  of 
-0.25  deg/day  as  the  rate  of  change  of  argument  of  perigee,  the  decay 
condition  was  then  expressed  as  the  intersection  of  two  curves  involving 
the  varying  value  of  perigee  and  the  varying  and  initial  values  of 
longitude  of  the  ascending  node. 

For  Molniya  satellites  in  orbits  of  63°  inclination,  [21]  found  that  the 
satellite  orbits  were  far  less  amenable  to  analytical  treatment,  because 
the  argument  of  perigee  could  librate  around  270°  and  the  inclination 
about  63.4°.  He  concluded  that  very  long  lifetimes,  on  the  order  of 
hundreds  of  years,  were  possible  if  the  argument  of  perigee  remained 

around  270°  for  cases  in  which  the  absolute  value  of  the  longitude  of  the 
ascending  node  was  below  60°.  In  1978  [22]  presented  a  graphical  method 


of  estimating  future  lifetimes  of  an  artificial  satellite  based  on  its 


current  rate  of  decay  and  by  using  an  exponential  atmosphere  model. 
Effects  of  the  asphericity  of  the  Earth,  effects  of  atmospheric  oblate¬ 
ness,  density  changes  due  to  the  11 -year  solar  cycle,  day-to-night,  and 
semiannual  variations  were  approximated  by  specified  correction  factors. 
Large  lunar-solar  perturbations  {greater  than  50  km)  were  handled  with 
numerical  integration.  For  high  eccentricity  orbits,  the  initial  value  of 
perigee  height  was  corrected  using  the  decay  rate  with  the  mass-to-area 
ratio  of  the  spacecraft.  These  methods  allowed  decay  date  predictions  up 
to  about  one  year  ahead  with  an  accuracy  of  about  10%  of  the  remaining 
lifetime . 

1.3.2  NORAD  Prediction  Models  and  the  Implementation  of  SDP4 

As  previously  mentioned,  NORAD  has  the  responsibility  of  maintaining  a 
catalog  of  element  sets  for  all  space  objects.  The  element  sets  are  then 
distributed  to  users  to  allow  prediction  of  the  future  position  and 
velocity  of  each  space  object.  The  satellite  observations  are  obtained  by 
the  Space  Detection  and  Tracking  System  (SPADATS)  network  of  "skin-track" 
radars  and  optical  sensors.  The  orbital  element  sets  are  based  on  statis¬ 
tical  processing  of  the  satellite  observations.  For  the  highly  eccentric 
orbits,  the  observations  are  primarily  obtained  from  the  deep  space  radars 
(Millstone  and  Altair)  and  from  the  GEODSS  optical  sensors. 

NORAD  has  utilized  several  propagation  models  for  prediction  of  satellite 


position  and  velocity.  They  will  be  summarized  here;  a  more  detailed 


description  is  available  in  [23].  The  earliest  model,  SGi,  was  developed 


in  1966  for  near  Earth  satellites.  It  simplified  the  work,  of  Kozai  for 


its  gravitational  model,  including  the  zonal  harmonics  J  >  and  J.  ,  and  it 


assumed  the  drag  effect  on  mean  motion  to  be  a  quadratic  in  time,  result¬ 


ing  in  a  cubic  variation  in  mean  anomaly  with  time.  The  drag  effect  on 


eccentricity  was  modeled  in  such  a  way  that  perigee  height  remained 


constant. 


A  second  propagation  model,  SGP4,  was  developed  in  1970  and  is  used  for 


near  Earth  satellites.  This  model  was  obtained  by  simplification  of  the 


more  extensive  analytical  theory  of  Lane  and  Cranford  which  used  the 


solution  of  Brouwer  for  its  gravitational  model  and  a  power  density 


function  for  its  atmospheric  model.  The  next  model  developed,  SDP4 ,  was 


an  extension  of  SGP4  to  be  used  for  deep-space  satellites  (which  NORAD 


defines  as  those  having  a  period  greater  than  225  minutes).  It  is 


described  in  more  detail  below. 


Another  propagation  model,  developed  but  never  implemented  operationally. 


was  SGP8.  The  SGP8  and  SDP8  models  have  the  same  gravitational  and 


atmospheric  models  as  the  SGP4  and  SDP4  models  respectively,  although  the 


form  of  the  solution  was  quite  different. 


Since  the  late  1970's,  SDP4  has  been  the  model  used  for  deep  space  satel¬ 


lites.  The  SDP4  theory  is  based  upon  the  restricted  four-body  solution 


for  resonating  satellites  without  drag  and  using  multiple  transformations 


developed  by  Hujsak  [24].  First,  the  four-body  oblate  Earth  problem, 
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transformations  are  int  reduced  an.i  the  transformed  dynami  m 1  sys'.'-r  is 
integrated  ana ly t i ca 1 1 y .  Me  then  used  the  method  of  averaqimj  t  r  derive  a 
multiply-transformed  dynamical  system  for  resonating  satellites  of  12-hr 
or  24-hr  orbital  periods.  This  dynamical  system  was  then  integrated 
analytically  except  for  a  numerical  evaluation  of  the  main  resonance 
effect. 


Hu  j  s  ak  used  the  spherical  harmonic  potential  of  the  Earth  in  terms  of 
associated  Legendre  functions  as  given  by  [25],  where  the  Earth  potential 
was  separated  into  the  zonal  potential  and  the  potential  due  to  sectoral 
and  tesseral  harmonics,  the  latter  of  which  was  expressed  in  terms  of  the 
F 1  mp  ( i  )  and  G  i  pq  ( e )  inclination  and  eccentricity  functions.  The 
specific  tesseral  resonance  terms  included  in  SDP4  for  the  12  hr  orbits 
are  of  degree  and  order: 

2.2  5,2 

3.2  5,4 
4,4 

The  gravitational  potential  of  the  Moon  or  Sun  was  represented  by  a  single 
term  of  the  expanded  series  of  Legendre  functions  of  the  phase  angle 
between  the  satellite  and  the  third  body. 

The  dynamical  system  was  a  function  of  two  fast  variables,  mean  anomaly 
and  Greenwich  sidereal  time,  and  satisfied  the  conditions  for  application 
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of  the  method  of  averaging.  The  resultant  first-order  short-periodic 
variations  due  to  Ji  are  given  in  Appendix  F  of  [24].  Second-order, 
short-periodic  variations  were  ignored.  Second  and  third  transformations 
were  then  applied  to  remove  the  mean  anomalies  of  the  Moon  and  Sun.  The 
triply-transformed  dynamical  system  for  non-resonating  satellites  was  then 
expressable  as  the  sum  of  zonal  expressions  and  third  body  effects  (see 
Appendices  D  and  E  of  [24])  which  varied  slowly  with  time. 

For  the  problem  of  resonance,  [24]  suggested  a  change  of  variables  which 
resulted  in  a  dynamical  system  which  was  a  function  of  a  single  fast 
variable,  namely,  the  sidereal  time.  The  singly-averaged  four-body  oblate 
Earth  potential  under  resonance  conditions  was  then  derived  using  the 
method  of  averaging.  A  resonating  triply-averaged  potential  was  then 
derived  after  averaging  over  the  lunar  and  solar  mean  anomalies. 

In  the  early  1990's  USAF  Space  Command  plans  to  convert  to  a  new,  more 
accurate  ephemeris  prediction  model  called  HANDE  [26]  to  predict  satellite 
position  and  velocity  and  to  maintain  the  National  Space  Surveillance 
Center  (NSSC)  catalog.  The  HANDE  model  was  developed  in  [27]  and  was 
designed  for  both  near-Earth  and  deep  space  satellites.  The  method 
includes  the  zonal  harmonics  ,  J3  ,  and  J4  for  its  gravitational  model 
and  a  dynamic  atmosphere  [28]  for  its  physical  atmospheric  model.  It 
includes  lunar  and  solar  effects  was  well  as  resonance  terms  for  satel¬ 
lites  with  1/2  day  and  1  day  periods.  HANDE,  as  all  the  NORAD  ephemeris 
models,  produces  predictions  in  an  Earth-centered  non-rotating  coordinate 


system  based  on  the  true  equator  and  mean  equinox  of  the  epoch  of  the 


element  set.  The  HANDE  theory  is  unique  relative  to  previous  NORAD 
theories  in  that  the  external  user  only  needs  the  ephemeris  reduction 
portion  of  the  theory. 

1.3.3  Semianalytical  Satellite  Theory 

Advanced  applications  of  space  for  both  military  and  scientific  missions 
require  fast,  reliable  orbit  predictions  for  artificial  satellites.  Orbit 
predictions  using  special  perturbation  theories,  based  upon  the  numerical 
integration  of  the  osculating  equations  of  motion,  are  accurate  but  slow, 
since  they  often  require  a  hundred  or  more  steps  per  satellite  orbit  to 
give  good  results.  Orbit  predictions  using  general  perturbation  theories 
are  fast  but  contain  large  errors  due  to  the  neglect  of  certain  perturba¬ 
tions.  The  SDP4  model,  for  example,  truncates  the  zonal  potential  after 
J4  ,  limits  the  third-body  potential  expansion  to  a  single  term,  assumes 
truncated  resonance  effects  based  on  a  12-hr  or  24-hr  orbit,  and  assumes  a 
power  density  function  for  its  atmospheric  model. 

Semiana ly tical  Satellite  Theory  (SST)  is  an  alternative  to  special  and 
general  perturbation  theories.  In  SST,  perturbations  that  can  be  express¬ 
ed  in  terms  of  a  disturbing  potential  have  that  potential  expressed  in 
nearly  singularity-free  equinoctial  elements.  The  perturbations  are  then 
put  in  Lagragian  Variation  of  Parameter  (VOP)  form.  This  allows  only 
small  deviations  from  Keplerian  motion  to  be  considered.  Those  perturba¬ 
tions  which  cannot  be  expressed  in  terms  of  potential,  such  as  drag  and 


solar  radiation  pressure,  are  expressed  in  Gaussian  VOP  equations 


SST  then  uses  the  Generalized  Method  of  Averging  [29,30]  to  separate  long- 


period  and  secular  terms  from  the  short-periodic  components  of  the  satel¬ 


lite  motion.  It  is  usually  sufficient  to  perform  this  separation  of  the 


so-called  "mean"  elements  from  the  short-periodic  terms  only  to  first 


order  in  small  parameters  of  the  perturbations;  however  some  second  order 


and  coupling  terms  have  been  included. 


The  analytical  or  numerical  averaging  technique  which  separates  these  mean 


components  results  in  a  smooth  mean  satellite  motion.  For  perturbations 


expressed  in  terms  of  a  potential  this  separation  is  primarily  accomplish¬ 


ed  by  averaging  all  the  perturbative  effects  over  a  fast  angular  variable 


for  a  satellite  period.  For  non-potential  perturbations  this  averaging 


takes  place  over  the  satellite  period  instead  of  a  cycle  of  an  angular 


variable. 


A  low-order  integration  technique,  e.g.,  the  fourth-order  Runge-Kutta 


method,  with  stepsizes  of  one  day  is  usually  sufficient  to  integrate  the 


averaged  satellite  motion  for  the  high  eccentricity  cases.  Interpolators 


are  used  to  produce  accurate  mean  orbital  elements  for  any  time  of  inter¬ 


est.  Reversing  the  averaging  transformation  results  in  the  short-periodic 


component  of  motion.  This  short-period  term  is  added  to  the  mean  motion 


to  obtain  the  satellite  osculating  state.  The  fact  that  the  calculation 


of  this  short-period  term  is  usually  only  required  for  the  particular  time 


of  the  desired  osculating  output  is  one  of  the  major  advantages  of  the 


SST.  This  is  because  the  mean  motion  already  allows  the  full  dynamics  to 


be  recovered.  The  short-period  portion  of  the  theory  also  incudes  an 
interpolator  concept  so  that  closely  spaced  output  requirements  can  be 
treated  efficiently. 

The  SST  used  in  this  investigation  has  been  implemented  in  the  CSDL 
version  of  the  Goddard  Trajectory  Determination  System  ( GTDS ) ,  a  multi¬ 
purpose  computer  system  originally  formulated  to  support  space  missions 
and  various  research  and  development  project  requirements  at  NASA  Goddard 
Space  Flight  Center  [31 ] .  The  Draper  version  includes  an  array  of  force 
models  that  can  be  called  upon  to  achieve  varying  levels  of  orbit  deter¬ 
mination  accuracy  for  different  satellite  orbits.  Shown  below  are  the 
force  models  developed  for  the  semianalytical  propagator. 

Mean  Element  Dynamics 

?  1 

Recursive  zonals  (closed  form  (c.f))  and  J2  e 
Recursive  tesseral  resonance  (en,  n>20) 

Recursive  lunar-solar  (single  averaging ) (c.f . ) 

Recursive  lunar-solar  (double  averaging ) (c.f . ) 

Solid  Earth  tide  (c.f.) 

Atmospheric  Drag  (and  J2 /drag  coupling ) (numerical ) 

Solar  radiation  pressure  (numerical) 

Short  Periodics 

2  1 

Recursive  zonals  (c.f.)  and  J2  e 
Recursive  tesseral  m-dailies  (c.f.) 

Recursive  tesseral  linear  combinations  (en,n>20) 
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Recursive  J2  secular/tesseral  m-dai ly  coupling  (c.f.) 

Recursive  lunar-solar  (c.f.) 

Atmospheric  drag  (numerical) 

Solar  radiation  pressure  (numerical) 

Numerous  studies  for  several  types  of  orbits  have  confirmed  the  accuracy 
of  the  theory  [32,33,34],  and  a  portable  version  has  also  been  developed 
[35]. 

1.4  THESIS  OVERVIEW 

The  primary  objective  of  this  thesis  is  to  test  the  semianalytical  satel¬ 
lite  theory  as  implemented  in  the  CSDL  version  of  the  GTDS  against  long 
arcs  of  real  data  of  highly  eccentric  orbits.  Of  primary  interest  will  be 
the  orbits  of  two  Soviet  Molniya  spacecraft  and  the  orbit  of  the  European 
Space  Agency  satellite  Exosat.  The  data  for  the  satellites  will  be  in  the 
form  of  NORAD  observations  and  element  sets.  Chapter  1  has  reviewed  the 
uses  of  elliptical  orbits  and  the  need  for  long-term  predictions.  It  has 
also  reviewed  the  application  of  artificial  satellite  theory  to  highly 
elliptical  orbits. 

Chapter  2  addresses  the  test  methodology  of  this  investigation.  The 
primary  topics  include: 

1.  Pre-processing  NORAD  data 

2.  Determining  an  initial  state  using  an  SST  DC 


3.  Generating  the  SST  prediction  file 

4.  Evaluating  the  SST  prediction 


Chapter  3  describes  the  results  obtained  using  real  data  for  specific 
orbits  using  the  methodology  of  Chapter  2.  The  data  available  for  this 
study  includes  observations  and  element  sets  for  the  following  objects: 


Object  Data  Type  Period 


NSSC  9829  (1 977-1 OA)  Mean  elements  February  1977  -  September  1986 

Molniya  2-17 


NSSC  14095  ( 1  983-51 A ) 
Exosat 


Mean  elements 
Observations 


June  1983  -  December  1985 
June  1983  -  December  1985 


NSSC  13964  (1983-25A)  Mean  elements 

Molniya  1-57  Observations 


May  1983  -  September  1985 
May  1983  -  December  1985 


The  rationale  for  the  dynamics  modelling  chosen  in  each  case  for  the 
differential  correction  and  for  the  semianalytic  propagator  is  given,  and 
an  analysis  of  the  comparison  and  difference  plots  is  made. 

Chapter  4  summarizes  the  results  of  the  previous  chapters  and  presents 


some  ideas  for  future  study 


CHAPTER  2 


TEST  METHODOLOGY 

The  primary  objective  of  this  thesis  is  to  test  the  semianaly tical  satel¬ 
lite  theory  as  implemented  in  the  CSDL  version  of  the  GTDS  against  long 
arcs  of  real  data  for  highly  eccentric  orbits.  Of  primary  interest  will 
be  the  orbits  of  two  Soviet  Molniya  spacecraft  and  the  orbit  of  the 
European  Space  Agency  satellite  Exosat.  The  data  for  the  satellites  is  in 
the  form  of  NORAD  element  sets  and  actual  observations.  Chapter  1  review¬ 
ed  how  elliptical  orbits  have  been  used  by  the  scientific  and  military 
communities,  described  the  need  for  long-term  predictions,  and  reviewed 
the  application  of  artificial  satellite  theory  to  highly  elliptical 
orbits.  Chapter  2  addresses  test  methodology. 

The  primary  topics  of  this  methodology  include: 

1 .  Pre-processing  NORAD  data 

2.  Determining  an  initial  state  using  an  SST  DC 

3.  Generating  the  SST  prediction  file 

4.  Evaluating  the  SST  prediction 

The  overall  framework  of  this  effort  is  shown  in  Figure  3a  and  Figure  3b. 
Figure  3a  illustrates  the  data  flow  in  which  a  relatively  short  arc  of 


sition  an 
v  data  in 


N  )RAP  element  sets  are  converted  to  mean  position  and  velocity.  The  mean 


position  and  velocity  data  are  used  as  lumped  observations  in  the  DC  to 
determine  initial  values  of  the  mean  elements  for  the  SST  prediction.  The 
results  of  this  prediction  are  compared  with  a  long  arc  of  NORAD  element 
sets  . 


An  alternative  approach  is  to  use  the  NORAD  observation  data  from  the  HDS 
to  determine  initial  values  for  the  SST  mean  element  set  (Figure  3b). 
Software  used  to  convert  range,  azimuth,  elevation,  and  range-rate  data 
from  NORAD  observation  format  to  GTDS  observation  format  had  been  built 
prevously.  For  this  study,  the  pre-processing  capability  was  expanded  to 
include  optical  data.  The  results  of  the  SST  prediction  are  .gain 
compared  with  a  long  arc  of  NORAD  element  sets. 

2.1  PRE-PROCESSING  NORAD  DATA  TO  GTDS  OBSERVATION  CARD  FORMAT 


Pr  jcessina  of  NORAD  data  proceeds  along  slightly  different  paths  depending 
on  whether  element,  sets  or  observations  are  used  in  the  DC.  Element  sets 
are  received  from  the  NORAD  Historical  Data  System  (HDS)  in  transmission 
carl  format  [23!  and  include  values  for: 
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XNODE 


Right  ascension  of  node  (degrees) 


EO 

Eccentricity 

OMEGA 

Argument  of 

perigee  (degrees) 

XMO 

Mean  anomaly 

(degrees  ) 

XNO 

Mean  motion 

(rev/day ) 

A  graphical  utility  written  in  IBM  FORTRAN  77  called  ADCEDIT  plots  each 
element  against  the  day  of  any  selected  year  of  data.  This  allows  "noisy" 
element  sets  to  be  editted. 

The  program  source  code  is  listed  in  Appendix  A.  The  routine  creates  a 
deck  of  data  cards  and  a  deck  of  control  cards  required  for  the  CSDL 
plotting  program  PLOT 4 B  [36].  Figure  4  is  an  example  that  plots  the 
argument  of  perigee,  mean  anomaly,  and  mean  motion  for  NSSC  9829  for 
1983.  The  element  set  that  caused  the  spike  was  consequently  editted  from 
further  pre-processing. 

The  next  step  in  the  pre-processing  of  element  set  data  is  to  generate 
position  and  velocity  in  GTDS  observation  format  [31 ] .  This  is 

accomplished  by  a  modified  version  of  the  NORAD  SDP4  routines  described  in 
[23].  DRIVER  has  been  replaced  by  the  utility  RUNADC,  and  subroutine  SDP4 
has  been  slightly  modified.  The  subroutines  ACTAN,  DEEP,  THETAG,  and 
FM0D2P  were  unchanged. 

The  algorithm  for  RUNADC  is  shown  in  Figure  5,  and  its  source  code  is 
listed  in  Appendix  A,  An  input  card  reads  the  satellite  designator,  the 
NORAD  generator  type,  and  the  data  card  format.  Error  messages  are  issued 


SV  9829 


NORAD  MEAN  ELEMENTS 


o  1983 


^3000  83040  83080  83120  83160  83200  83240  83280  83320  63360 

DATE  ( YYDDD  1 


H33000  83040  83080  83120  83160  83200  83240  83280  83320  83360  63400 

DATE  I YYDDD) 


B300C  83040  83080  83120  83160  83200  83240  83280  83320  63360  83400 

DATE  (YYDDD) 

Figure  4.  NSSC  9829  NORAD  Element  Sets-1983  (OMEGA,  XMO ,  XNO) 
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Figure  5.  Subroutine  RUNADC  Algol  ithm  for  Pre-Proeessi ng  NORAI  Element  Set 


if  the  generator  type  is  not  'SDP4*  or  the  data  card  format  is  other  than 


'TRANS'.  The  input  card  also  contains  the  type  of  output  desired  (mean  or 
osculating),  the  start  date  and  end  date  for  the  GTDS  observation  card 
file  and  another  temporary  dataset  of  NORAD  elements,  a  search  interval, 
and  a  suppression  interval. 

The  first  NORAD  element  card  is  read  and  examined  to  see  if  its  epoch 
occurs  after  the  start  date  specified  on  the  input  card.  When  a  satisfac¬ 
tory  epoch  is  found,  the  data  is  examined  to  see  if  the  next  epoch  occurs 
within  the  suppression  interval,  typically  12-15  hours,  specified  on  the 
input  card.  The  suppression  interval  helps  ensure  that  the  element  set  is 
the  final  element  set  in  a  NORAD  differential  correction.  The  element  set 
is  then  input  into  the  routine  SDP4.  "TSINCE"  is  set  equal  to  zero 
because  the  SDP4  subroutine  will  not  be  used  to  generate  predictions. 
That  is,  the  position  and  velocity  observation  data  are  generated  by  SDP4 
at  the  time  epoch  attached  to  the  input  element  set.  Also,  the  J2  short- 
period:"  model  included  in  the  original  SDP4  may  be  circumvented  by 
specifying  'mean'  as  the  type  of  output  motion  desired.  The  utility  thus 
generates  single-averaged,  or  'mean'  data.  This  simplifies  the  GTDS  DC 
process  because  only  the  mean  dynamics  need  to  be  included.  These  mean 
position  and  velocity  data  are  written  in  GTDS  observation  card  format  and 
use',  to  calculate  Keplerian  and  equinoctial  elements.  These  data  are 
written  into  a  temporary  data  set  of  "NORAD"  points  which  will  t>e  later 
plot  r«>d  with  tne  "  GTDG "  predictions.  The  next  element  set  accepted  by  th« 


pre-processor  will  occur  after  the  search  interval,  typically  5-7  days. 
This  search  interval  was  necessary  to  limit  the  number  of  GTDS  observation 
cards.  Element  sets  are  accepted  until  the  specified  end  date  is  reached. 

The  observation  cards  thus  created  serve  as  the  input  to  the  GTDS  differ¬ 
ential  correction  program  which  is  used  to  determine  the  initial 
conditions  for  the  ephemeris  prediction  using  the  semianalytical  theory. 

The  procedure  of  using  the  element  sets  to  initialize  the  DC  has  the 
advantage  that  the  position  and  velocity  observations  require  little 
auxiliary  data  (such  as  observation  noise  and  bias,  station  location, 
etc.).  The  disadvantage  of  this  procedure  is  that  the  elements  sets 
contain  errors  originating  in  the  original  observations  and  in  the 
processing  used  to  generate  the  element  sets.  Determination  of  the 
initial  SST  mean  elements  may  require  the  NORAD  element  set  error 
statistics,  and  this  data  is  not  maintained  in  the  HDS . 

When  actual  NORAD  observations  are  used  to  determine  the  initial  condi¬ 
tions  for  the  DC,  pre-processing  is  necessary  to  transform  the  data  from 
NORAD  to  GTDS  format.  For  radar  observations,  which  may  include  range, 
azimuth,  elevation,  and  range  rate,  the  pre-processing  is  primarily  a 
matter  of  converting  time  and  units.  For  the  optical  data,  the  NORAD  HDS 
provides  the  right  ascension  and  declination  data  in  a  true  equator  mean 
equinox  of  date  coordinate  system.  This  data  must  be  transformed  to  the 
Mean  of  195'i.e  coordinate  system  before  it  can  be  used  with  the  GTDS 


Pre-processing  of  the  observations  is  performed  by  the  CSDL  utility 
RUNADCOB.  For  optical  observations,  RUNADCOB  utilizes  another  subroutine 
ASTRON  to  convert  to  the  Mean  of  1950.0  coordinate  system.  The  precession 
and  nutation  matrices  for  the  coordinate  transformation  are  calculated  by 
the  subroutine  PRENUT  whose  algorithms  originated  in  the  NORAD  SACEPH 
program.  Listings  of  RUNADCOB,  ASTRON,  and  PRENOT  are  found  in  Appendix 
A.  The  observation  cards  produced  by  RUNADCOB  serve  as  the  input  to  the 
GTDS  DC  program. 

2.2  SST  DIFFERENTIAL  CORRECTION  PROCESS 

The  primary  purpose  of  the  GTDS  Differential  Correction  (DC)  program  is  to 
estimate  the  satellite  orbit  and  auxilliary  parameters  (such  as  the  drag 
coefficient,  solar  pressure  coefficient,  etc.).  The  estimation  algorithm 
used  in  the  DC  program  is  called  the  weighted  lest  squares  with  a  priori 
algorithm  or  the  Bayesian  weighted  least  squares  algorithm  [31].  It 
minimizes  the  sum  of  the  squares  of  the  weighted  residuals  between  actual 
and  computed  observations,  while  simultaneously  constraining  the  model 
parameters  to  satisfy  the  a  priori  conditions  to  within  a  specified 
uncertainty.  Both  first-  and  second-order  statistics  (i.e.,  the  mean  and 
covariance  matrices)  are  determined  for  the  estimated  variables. 

The  user  provides  the  DC  program  with  the  initial  estimate  of  the  solve- 
for  vector  at  a  specified  epoch.  (For  this  effort,  the  initial  estimate 


for  the  elements  was  taken  from  the  NORAD  element  set  data.)  The  program 


can  accept  a  variety  of  observation  types  in  one  of  several  input 


coordinate  systems.  Observation  uncertainties,  atmospheric  model,  and 
spacecraft  area  and  mass  parameters  may  also  be  specified.  The  DC  program 
can  access  the  complete  semianalytical  theory  as  described  in  Section 

1.3.3. 


2.3  SST  EPHEMERIS  PREDICTION  PROCESS 

The  function  of  the  ephemeris  generation  program  GTDS  EPHEM  is  to  compute, 
from  prescribed  initial  conditions,  the  value  at  a  specified  time  of  the 
vehicle  state  and,  optionally,  the  state  partial  derivatives.  In  order  to 
meet  varying  precision  and  efficiency  requirements,  several  orbital 
theories  have  been  provided,  ranging  from  a  first-order  analytic  theory  to 
a  high-precision  Cowell-type  numerical  integration.  The  state  partial 
derivatives  can  be  computed  by  precision  numerical  integration  of  varia¬ 
tional  equations.  The  state  partial  derivatives  with  respect  to  the 
initial  state,  i.e.,  the  state  transition  matrix,  can  optionally  be 
generated  by  a  two-body  analytic  approximation. 

GTDS  EPHEM  includes  the  complete  semianalytical  theory  as  described  in 
Section  1.3.3.  The  user  selects  the  appropriate  model  for  the  orbit  under 
study  and  inputs  the  initial  elements  at  epoch.  The  GTDS  EPHEM  program 
can  accept  the  input  conditions  in  one  of  several  coordinate  systems.  A 
desired  atmospheric  model,  spacecraft  area,  and  mass  parameters  may  also 
be  specified. 


The  GTDS  EPHEM  program  writes  the  prediction  file  of  mean  position  and 


velocity  data  into  an  0RB1  file.  The  FORTRAN  utility  RD0RB1  reads 
position  and  velocity,  calculates  the  corresponding  Keplerian  and 
equinoctial  elements,  and  writes  this  "GTDS"  data  into  a  file  with  the 
same  format  as  the  file  of  "NORAD"  points.  A  copy  of  RD0RB1  is  included 
in  Appendix  A. 


2.4  PLOTTING  ELEMENT  COMPARISONS  AND  DIFFERENCES 

Comparison  and  difference  plots  of  the  "GTDS"  predictions  and  the  "NORAD" 
points  are  created  using  the  CSDL  plotting  program  PLOT 4 B  [36].  Because 
of  limitations  of  PL0T4B  it  was  necessary  to  interpolate  the  "GTDS"  pre¬ 
dictions  (spaced  regularly  in  time)  to  the  "NORAD"  output  from  the 
modified  SDP4  routine  (spaced  irregularly  in  time).  This  interpolation  is 
performed  by  program  PLOTTER. 

The  FORTRAN  utility  PLOTTER  reads  the  set  of  position  and  velocity 
components,  the  Keplerian  elements,  and  the  equinoctial  elements  from  the 
temporary  dataset  of  "GTDS"  predictions  created  by  the  RD0RB1  utility. 
It  first  reads  the  input  control  cards  containing  the  satellite  designa¬ 
tor,  the  time  interval  at  which  the  GTDS  predictions  are  spaced,  the 
comparison  start  date,  the  start  day  of  the  plot  and  the  final  day  of  the 
plot  in  number  of  days  since  the  start  date,  the  FORTRAN  file  numbers  of 
the  two  data  sets,  and  the  gravitational  constant. 
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The  angular  velocities  of  mean  anomaly  and  mean  longitude  are  used  to 
convert  the  "sawtooth"  values  of  longitude  of  ascending  node,  argument  of 
perigee,  mean  anomaly,  and  mean  longitude  to  linearized  values.  These 
linearized  values  for  the  angular  elements  can  then  be  used  to  perform  a 
five-point  Lagragian  interpolation  of  the  SST  predictions  to  the  times 
corresponding  to  the  NORAD  mean  elements.  The  numerical  difference 
between  each  SST  prediction  and  NORAD  mean  element  is  calculated,  and 
first-order  statistics  (mean  difference  and  standard  deviation)  are 
computed  for  the  entire  comparison  interval. 

The  radial,  cross-track,  and  along-track  errors  using  the  "NORAD"  points 
as  the  estimated  trajectory  and  the  "GTDS"  predictions  as  the  true 
trajectory  are  then  calculated.  Unit  vectors  for  these  errors  are  defined 
by: 

Radial  error  H 

Cross-track  error  C 

Along-track  error  L 

The  program  then  writes  PL0T4B  control  cards  and  data  cards  for  comparison 
plots,  difference  plots,  and  error  plots.  A  copy  of  the  source  code  for 
PLOTTER  is  found  in  Appendix  A. 
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CHAPTER  3 


TEST  RESULTS 


The  primary  objective  of  this  thesis  was  to  test  the  semianalytical  satel¬ 
lite  theory  as  implemented  in  the  CSDL  version  of  the  GTDS  against  long 
arcs  of  real  data  for  highly  eccentric  orbits.  Chapter  2  addressed  the 
test  methodology  used  in  this  investigation  and  described  the  software 
development  that  was  undertaken  to  interface  the  NORAD  data  to  the  GTDS 
format  and  then  to  interface  GTDS  predictions  and  NORAD  elements  to  a 
simple  plotting  program  package.  The  present  chapter  summarizes  the 
comparisons  between  the  SST  predictions  and  the  NORAD  real  data  using  a 
9-year  arc  of  element  sets  for  NSSC  9829,  the  Soviet  Molniya  2-17  space¬ 
craft,  an  18-month  arc  of  element  sets  for  NSSC  14095,  the  ESA  Exosat,  and 
an  18-month  arc  of  observations  for  NSSC  13964,  the  Soviet  Molniya  1-57 
spacecraft.  Comparison  and  difference  plots  are  given  for  all  three 
satellites . 


3.1  INITIALIZATION  OF  SST  USING  ELEMENT  SETS 


3.1.1  NSSC  9829 


3. 1.1.1  Mission  and  Operations 

The  NSSC  9829  spacecraft,  Molniya  2-17,  was  launched  11  February  1977. 
The  spacecraft  entered  a  supra-semi -synchronous  transfer  orbit  with  a 
period  of  735  minutes  [1]  and  was  allowed  to  drift  westward  to  its 
operational  location.  Its  initial  operational  orbital  parameters  on  1  May 
1977  [1]  are  shown  in  Table  1. 


Table  1 .  NSSC  9829  Initial  Operational 
Orbit  Elements 

Element 

Value 

Semimajor  axis 

Apogee  height 

Perigee  height 

Eccentricity 

Inclination 

Argument  of  perigee 
Period 

26554  km 

39853  km 

498  km 

0.741 0 

62.9° 

280° 

717.67  minutes 

NORAD  Coordinates 

The  analysis  for  this  satellite  was  based  solely  on  element  set  data. 
However,  a  review  of  the  graphs  created  by  ADCEDIT  indicated  there  were  no 
apparent  maneuvers  during  the  period  of  the  available  NORAD  data. 

Figure  4  is  a  plot  of  argument  of  perigee,  mean  anomaly,  and  mean  motion 
for  NSSC  9829  for  1983.  The  "spikes"  on  all  three  plots  clearly  indicate 
a  bad  element  set.  These  "spikes”  typify  the  element  sets  that  were 
deleted  from  further  pre-processing  and  comparisons.  These  element  sets 
are  listed  m  Tabl»  2. 


Table  2.  NSSC  9829  Element  Set  Edits 


Element  set  date 

Questionable  Element 

( YYDDD .DDDDDDDD ) 

77191 .3281601 3 

XNDT2 

77196.31 168989 

XNDT2 

77298.97607433 

BSTAR 

77312.93081504 

XNDD6 

7921 1 .99956568 

XNDT2 ,  XINCL 

79271 .81988559 

XNDT2 

81256.60694402 

XINCL 

81  273.5553391  3 

BSTAR 

81275.54929572 

BSTAR 

82228.57192200 

XNDT2,  XINCL,  BO,  OMEGA 

82245.01890805 

XNDT2,  BSTAR 

83259.82751 893 

XNDT2,  XINCL,  EO,  OMEGA,  XMO,  XNO 

It  should  be  noted  that  the  percentage  of  element  sets  editted  was 
extremely  small  (about  1%  of  the  9  year  arc). 


3. 1.1. 2  SST  DC  and  Prediction 


The  initial  estimate  for  the  vehicle  state  was  taken  from  a  NORAD  mean 
element  set.  Six  months  of  data  were  used  in  the  fit,  and  the  epoch  was 
chosen  to  minimize  noise  in  the  mean  motion  rate  and  to  ensure  perigee  was 
outside  the  main  portion  of  the  atmosphere.  Table  3  summarizes  the  force 
models  employed  in  the  SST  DC  and  prediction.  The  initialization  proce¬ 
dure  automatically  computes  the  maximum  powers  of  a/R  and  e  to  be  included 
in  the  power  series  expansions  [37],  Drag  terms  were  not  included  because 
perigee  height  was  outside  the  main  portion  of  the  atmosphere,  and  solar 
pressure  terms  were  not  included  because  early  rsns  showed  their 
contribution  was  minimal. 


Table  3.  Force  Model  used  for  NSSC  9829 


Zonals:  J2  thru  Jj  g  ,  e  ,  GEM  9 

Tesseral  resonance:  GEM  9,  (2,2)  through  (6,6),  e3 
Non-central  bodies: 


Moon  (a/r)9,  e7 


Sun  (a/r)3,  e3 


Solve-for  parameters:  mean  equinoctial  elements 
Drag:  off 

Solar  pressure:  off 


The  a  priori  values  and  the  final  values  for  the  Keplerian  elements  and 
their  standard  deviations  are  shown  in  Table  4.  The  assumed  observations 
standard  deviations  were  20  km  for  each  position  coordinate  and  20  m/sec 
for  each  velocity  coordinate. 


Table  4.  Results  of  DC  tor  NSSC  9829 


Epoch 

Aug  3,  1979, 

23  hours,  42  min 

12.000  sec 

Element 

A  priori  Value 

Final  Value 

Standard 

Deviation 

a  (km ) 

26556.8740 

26556.9582 

1 .67468E-03 

e 

0.6992405 

0.6990986 

8.42488E-05 

i  ( 0  ) 

63.1 44666 

63.173001 

0.59929 

n  (°) 

190.635924 

190.619681 

0.42451 

u  (°) 

281 .54328 

281 .59624 

1.2282 

M  (  0  ) 

13.30040 

1 3.2931 5 

0.201 59 

We  1  gh  ted 

RMS:  0.2893 

Mean  of  1950.0  Coordinates 
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The  final 


values  of  the  Keplerian  elements  and  the  epoch  indicated  in 


1 


Table  4  were  used  as  initial  conditions  for  an  S3T  prediction.  The  same 
force  models  shown  in  Table  3  were  used  to  create  an  ephemeris  prediction 
file  of  about  2330  days. 

3. 1.1. 3  Comparison  and  Difference  Plots 

The  results  of  comparing  of  the  "NORAD*  points  to  the  "GTDS"  predictions 
using  the  GEM  9  Earth  gravitational  model  and  a  6  x  6  harmonic  field  are 
shown  in  Figures  6-20.  Figures  6-12  are  comparison  plots  of  the  "NORAD" 
and  "GTDS"  predictions  for  the  Keplerian  elements  for  the  approximately 
2330  days  of  available  NORAD  elements  following  the  indicated  epoch  of 
August  1979.  The  mean  difference  and  standard  deviation  after  326  com¬ 
parisons  for  the  Keplerian  and  equinoctial  elements  are  shown  in  Table  5. 


Table  5.  Results  of  C 

Comparisons  for  NSSC 

9829 

Element 

Mean  Difference 

Standard 

Deviation 

Semimajor  axis  (km) 

.25046 

.19496 

Eccentricity 

0.22437E-03 

0 . 38909E-03 

Inclination  ( 0  ) 

0. 33941 E-01 

0 . 87707E-01 

Longitude  of  ascending  node  (°) 

0.22964 

0.13148 

Argument  of  perigee  (°) 

0.11 699 

0.10242 

Mean  anomaly  (°) 

1 .0960 

0.78569 

Radius  of  perigee  (km) 

5.9652 

10.348 

Radius  of  apogee  (km) 

5.9333 

10.242 

H 

0.961 50E-03 

0. 77753E-03 

K 

0.98751 E-03 

0.70091 E-03 

P 

0. 1 6495E-02 

0.1 7350E-02 

2 

0. 1 6255E-02 

0.1  1 163E-02 

Mean  longitude  (°) 

0.99281 

0.74933 

Norad  Coordinates 

VVWVT  Fry*  yrji  W'/TT^ 


The  slowly  varying  Keplerian  element  histories  (Figures  6-10)  and  the  mean 
anomaly  difference  (Figure  18)  clearly  demonstrate  the  ability  of  the 
semiana lytical  theory  to  predict  the  dominant  motion  experienced  by  the 
NORAD  element  sets. 


Gravitational  models  and  tesseral  resonance  terms  other  than  those  listed 
in  Table  3  were  used  in  various  DC  runs.  In  each  case  the  semimajor  axis 
difference  exhibited  a  frequency  very  similar  to  the  resonance  period 
exhibited  by  the  semimajor  axis  difference  plot  (Figure  13).  A  marked 
improvement  was  noted  when  a  6x6  field  with  GEJ4  9  coefficients  was  used  to 
generate  the  semianalytical  prediction  in  Figures  6-12. 


When  an  8x8  field  with  GEM  9  coefficients  was  used,  the  DC  converged  very 
unsatisfactorily  with  an  overall  weighted  RMS  value  of  2.071  over  the  6 
month  fit  span  (comi  -re  this  to  the  weighted  RMS  value  of  0.2893  shown  in 
Table  4). 


Consideration  of  this  poor  fit  led  to  several  Precise  Conversion  of 
Elements  (PCE)  runs  by  R.  Proulx  (CSDL)  which  generated  predictions  using 
a  Cowell  numerical  integration  scheme.  These  predictions  very  closely 
resembled  the  SST  predictions  using  the  8x8  field.  It  was  at  this  point 
that  usage  of  actual  NORAD  observations  to  improve  the  fit  was  considered. 
This  led  to  the  method  of  initialization  of  the  SST  prediction  using 
observation  data  which  is  described  in  Section  3.2. 
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Table  6.  NSSC  14095  Initial  Operational 
Orbit  Elements 


Element 


Value 


Apogee  height 
Perigee  height 
Inc  1 1 na t i on 
Argument  of  perigee 
Per  i  ocl 


2  X  10’  km 
500  km 
72.5° 

28b. 5° 

99  hours 
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Table  7.  NSSC  14095  Element  Set  Edits 

Element  set  date 
( YYDDD.DDDDDDDD ) 

Questionable  element 

83276.9999999 

85167.9999999 

85293.9999999 

XNDT2,  XNDD6,  XINCL,  XNODE,  EO,  OMEGA 
XNDT2,  XNDD6,  XINCL,  XNODE,  EO,  OMEGA 
XNDT2,  XNDD6,  XNO 

3. 1.2. 2  SST  DC  and  Prediction 

As  was  the  case  for  NSS<'  9829,  the  initial  estimate  for  the  vehicle  state 
for  NSSC  14099  was  taken  from  a  NORAP  mean  element  set.  Six  months  of 

data  wer"  used  in  the  fit,  and  epoch  was  chosen  to  ensure  a  penqee  out¬ 
side  tiit-  main  (x:>r*ion  of  Uie  atmosphere.  Table  8  lists  the  force  models 
uM  li.'.ei  i  ri  the  SST  DC  and  prediction.  The  maximum  powers  of  a  'V  and  e 

ire  lmled  i  r .  the  j»»wet  series  expansions  are  larger  than  those  included  for 

Ns  ■  '*•<.’  *  i  Ta!  '.*•  1'  he  a  use  of  *  he  limner  period  and  higher  ecrontritity 

•  ti  »  hx"sa*  hr. in  effects  were  not  included  because  penqee 

he.,]',*  wa-  >  -i»  ••  i  de  the  i f  mosphe  r  e ,  and  solar  pressure  effects  were  shown 
•  be  mi  "  :  -si  ;  . 

Table  8.  Pore,.  Model  used  for  NSSC  14h9‘. 

Ilona  Is:  d  .  through  e^,  GEM  9 

Tessera  1  r  e  s  ona  rnr  e  :  off 

N- >r  — -ent  r  a  1  bodies: 

Mi " >n  !  a  'r  ) *  1  ,  e 1  M 

Sun  ' a  r  ) ^ ,  e^ 

.1  ve  -f  u  parameter*.:  mean  ixpiinoe*  ial  elements 


The  a  priori  values  and  the  final  values  for  the  Keplerian  elements  and 
their  standard  deviations  are  shown  in  Table  9.  The  assumed  observation 
standard  deviations  were  40  km  for  each  position  coordinate  and  40  m/sec 
for  each  velocity  coordinate. 


Table  9.  Results  of  DC  for  NSSC 

14095 

Epoch : 

Aug  19,  1983,  23  hours,  59  min 

59.000  sec 

Element 

T  priori  Value 

Final  Value 

Standard 

Deviation 

a  (km  ) 

102351  .2699 

102368.2782 

0.25212 

e 

0.928471 

0.9260O1 

1  .5261 6L-4 

i  ( 0 ) 

71  .98431 

72.12634 

1 .2160 

•'  ( 0 ) 

185.5514 

184.9844 

4.3266 

(  0  ) 

283.5979 

283.6980 

1 . 4847 

M(  o  ) 

233.4717 

233.641 1 

2.6055 

We  l  qh  ted 

RMS:  2.357 

Mean  of  1990.0  Coordinates 

It  is  interesting  to  compare  the  standard  deviation  obtained  for  the  mean 
semi  roa  'or  axis  ( 2'-2  meters)  with  the  standard  deviation  obtained  f  or  tin1 
mean  semima’or  axis  f  ■  >r  NS  SC  98  (  1  .».?  meters). 

The  final  values  of  the  Keplerian  elements  and  the  epoch  indicated  i  ;i 


Table 

4  were 

Used 

a  s  initial 

cond l t l ons 

for  an  SST  prediction. 

The  same 

fere*.' 

mode  i  s 

shown 

m  Table  8 

were  used 

to  create  an  ephemeris 

predict : on 

f  1  1  e 

of  abf>u' 

■  H4 

da  vs  . 

3. 1.2. 3  Comparison  and  Difference  Plots 


Figures  23-29  are  comparison  plots  of  the  "NORAD"  and  "GTDS"  predictions 
for  the  Keplerian  elements  for  the  approximately  250  days  of  available 
NORAD  data  following  the  indicated  epoch  of  August  1983. 


This  prediction  length  was  chosen  because  plots  of  the  longer  prediction 
length  (840  days)  revealed  two  distinct  changes  in  slope  of  the  comparison 
plot  for  mean  anomaly  (Figure  22).  This  suggested  that  the  Exosat 
performed  maneuvers  during  May  1984  and  February  1985.  The  mean 

difference  and  standard  deviation  after  27  comparisons  for  the  Keplerian 
and  equinoctial  elements  are: 


Table  10.  Results  of 

Comparisons  for  NSSC  14095 

Element 

Mean  Difference 

Standard 

Devi ation 

Semimaior  axis  (km) 

77.274 

1 49.98 

Eccentr .cl ty 

0.1 2844E-02 

0.73346E-C3 

Inc  1 1  nat.  i  on  (  0  ) 

0.10235 

0.75591 E-01 

Longitude  of  ascendinq  node  (°) 

0.  38547 

0.27984 

Argument  of  perigee  (°) 

0.12361 

0.10469 

Mean  anomaly  (°) 

0.28861 

0.53563 

Radius  of  perigee  (knl 

133.69 

82.289 

Radius  of  apogee  (km) 

185.80 

248.62 

H 

0. 1  2355E-02 

0. 1 0365E-02 

K 

0 . 42946E-02 

0. 32847E-02 

P 

0 . 46902E-02 

0. 33793E-02 

c 

0.1 7882E-02 

0.1 4650E-02 

Mean  longitude  (°) 

0.46272 

0.47961 

Figures  30-37  are  the  plots  of  the  difference  between  the  "GVDS" 
predictions  and  the  "NORAD"  points  for  the  Keplerian  elements  plotted 
against  the  number  of  days  since  August  1983. 
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Figure  29.  NSSC  14095  Radius  of  Apogee  Comparison 
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3. 2. 1.1  Mission  and  operations 

The  NSSC  13964  spacecraft,  Molniya 
replacement  fer  Molniya  1-52.  Its 
700  minutes  [40],  indicating  the 
operational  location.  Its  initial 
April  1983  are  shown  in  Table  11. 


1-57,  wart  launched  2  April  1983  as  a 
initial  transfer  orbit  had  a  period  of 
spacecraft  drifted  eastward  to  its 
operational  orbital  elements  [9]  for  7 
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Table  11 


NSSC  13964  Initial  Operational 
Orbit  Elements 


Element 

Value 

Apogee  height 

Perigee  height 
Inclination 

Period 

39,877  km 

482  km 

62.93  degrees 
717.85  minutes 

3. 2. 1.2  SST  DC  and  Prediction 

The  NORAD  observation  data  was  comprised  of  radar  observations  and  optical 
data.  The  observation  weights  that  were  used  are  given  in  Tables  12  and 
!  13.  These  result  from  [41  ,  42,  43],  Also,  the  range  and  range-rate 

weights  have  been  heuristically  adjusted  to  consider  the  particular 
short-periodic  model  chosen  for  use  in  these  semianalytical  DC  runs.  A 
NORAD  mean  element  set  was  used  as  an  initial  estimate  for  the  DC  run. 
Processing  included  the  first  three  months  of  1985  observations.  Tables 
14  and  15  list  the  force  models  utilized  in  the  SST  DC  and  prediction.  It 
was  necessary  to  include  short-periodics  modelling  because  the  real  data 
was  in  the  form  of  actual  observations. 


Table  12.  Radar  Observatic 

jn  Data  Sigmas 

! 

Station 

Range 

Azimuth 

Elevation 

Range -rate 

(m ) 

(mdeg ) 

(mdeg ) 

(cm/s  ) 

Eglin 

67. 

17. 

15. 

n/a 

(399) 

Altair 

60. 

31  . 

17. 

15. 

(666) 

Mi 11s  tone 
(369) 
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Table  13.  Optical  Observation  Data  Sigmas 


Station 

Right  Ascension 
(mdeg ) 

Declination 
(mdeg ) 

Primary  Telescope 
(231 , 

232, 

221 , 

222, 

211 , 

212) 

5.5 

5.5 

Secondary  Telescope 
(233, 

223, 

213) 

6.9 

6.9 

Table  14.  Mean  Dynamics  Model  used  for  NSSC  13964 
Zonals:  J2  through  Jg  ,  e5  ,  WGS72  (12x12) 

Tesseral  resonances  WGS72  (12x12),  (2,2)  through  (8,8),  e20 
Non-central  bodies: 

Moon:  (a/r)9,  e7 

Sun:  (a/r  )3  ,  e3 

Solve-for  parameters:  mean  equinoctial  elements, 

solar  radiation  coefficient 

Drag:  off 

Solar  pressure:  on  (Area  =  3  rr£  ,  mass  =  250  kg) 


Table  15.  Short-Periodics  Model  used  for  NSSC  13964 


Zonals:  J2  through  ,  WGS72  (12x12) 


Tesseral  m-dailies:  4x4 


J2  secular /m -daily  coupling:  off 
High-frequency  tesserals:  off 


Non-central  bodies: 
Moon:  (a/r  J1  0  ,  e1  0 

Sun:  (a/r)4,  e4 


Solar  pressure:  off 
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Figure  40.  NSSC  13964  DC  Residual  (Elevation) 
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Figure  41.  NSSC  13964  DC  Residual  (Range  Rate) 
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Figure  42.  NSSC  13964  DC  Residual  (Right  Ascension) 
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Epoch 


Element 


A  priori  Value 

Final  Value 

26554.300 

26554.302 

0.7179000 

0.7179169 

63.09000 

63.09183 

163.31300 

163.31333 

279.65000 

279.65431 

12.627000 

12.627693 

1  .2000000 

1 .3425221 

4. 371 95E-05 
2.61 725E-07 
1  .02427E-02 
1.2411 3E-02 
7 . 33589E-03 
2. 28709E-03 
0.758E-03 


Residual  Plot  Symbol 
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mm 


Table  18.  Observation  Summary  by  Type  for  NSSC  13964 


Type 

Total  Number 

Number  Accepted 

Range 

1  31 

73  (55%) 

Azimuth 

131 

100  (76%) 

Elevation 

1  31 

107  (81%) 

Range  Rate 

102 

68  (66%) 

Right  Ascension 

501 

434  (86%) 

Declination 

501 

474  (94%) 

A  comparison  of  Table  16  and  Table  4  demonstrates  how  the  use  of  observa¬ 
tion  data  improved  the  determination  of  initial  conditions.  The  standard 
deviation  for  the  mean  semimajor  axis  improved  from  a  value  of  1.67m  to  a 
value  of  0.0437m. 

3. 2. 1.3  Comparison  and  Difference  Plots 

The  final  values  of  the  Keplerian  elements  and  the  solar  radiation  coeffi¬ 
cient  indicated  in  Table  16  were  used  as  initial  conditions  for  an  SST 
prediction.  The  force  model  of  Table  14  was  used  to  create  an  ephemeris 
prediction  file  of  about  250  days.  Figures  44-50  are  comparison  plots  of 
the  "NORAD"  and  "GTDS"  predictions  for  the  Keplerian  elements.  The  mean 
differences  and  standard  deviations  after  126  comparisons  are  shown  in 
Table  19. 

Figures  51-58  are  the  plots  of  the  difference  between  the  "GTDS"  predic¬ 
tions  and  the  "NORAD"  points  for  the  Keplerian  elements.  Note  that  the 
comparison  difference  plot  for  semimajor  axis  (Figure  51 )  exhibits  none  of 
the  resonance  period  error  exhibited  by  NSSC  9829  using  element  sets 


© 


7 


(Figure  13).  The  remaining  difference  plots  (Figures  52-58)  show  no 
secular  trends,  with  the  exception  of  argument  of  perigee  (Figure  55). 


Table  19.  Results  of 

Comparisons  for  NSSC  13964 

Element 

Mean  Difference 

Standard 

Deviation 

Semimajor  axis  (m) 

67.695 

64.706 

Eccentricity 

0. 76003E-04 

0.86377E-04 

Inclination  (°) 

0.78965E-02 

0.82670E-02 

Longitude  of  ascending  node  (°) 

0.1 5868E-01 

0.1 0429E-01 

Argument  of  perigee  (°) 

0.1 7679E-01 

0.11 293E-01 

Mean  anomaly  (°) 

0.19882 

0.1 3963 

Radius  of  perigee  (km) 

2.0123 

2.2880 

Radius  of  apogee  (km) 

2.061 1 

2.3298 

H 

0.1 1456E-03 

0.11 290E-03 

K 

0.1 9659E-03 

0.1 2057E-03 

P 

0.1 7967E-03 

0.11 745E-03 

<2 

0.89706E-04 

0.781 50E-04 

Mean  longitude  (°) 

0.19429 

0.1 3477 

The  magnitudes  of  the  differences  between  the  SST  predictions  and  the 
NORAD  points  are  much  less  for  the  case  where  observation  data  were  used. 
(Compare  Table  19  with  Table  5.)  This  clearly  demonstrates  that  the 
ability  of  the  semianalytical  theory  to  predict  dominant  motion  is 
enhanced  by  initialization  with  observation  data. 


3.2.2  Simulation  of  SDP4  with  SST  for  NSSC  13964 


An  effort  was  made  to  account  for  the  errors  between  the  SST  predictions 
and  the  NORAD  mean  elements  seen  for  the  Molniya  orbits.  This  was 
accomplished  by  configuring  the  SST  with  the  dynamic  modelling  used  in 
SDP4  and  performing  a  DC  over  a  short  interval  of  observation  data  for 
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Figure  54.  NSSC  13964  Longitude  of  Ascending  Node  Difference 
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Figure  55.  NSSC  13964  Argument  of  Perigee  Difference 
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Table  20.  Mean  Force  Models  in  SDP4  Simulation 


Zonals:  through  J4  ,  e2  ,  WGS  72  (12x12) 

Tesseral  Resonance:  WGS  72  (12x12),  (2,2),  (3,2),  (5,2), 

(4,4),  (5,4),  e20 

Non-Central  Bodies: 

Moon:  (a/r)2,  e2 
Sun:  (a/r)3,  e3 

Drag:  off 

Solar  Pressure:  off 


Table  21 . 

Short-Periodic  Force  Models  in  SDP4  Simulation 

Zonals:  J2 

only 

M-daily  Tesserals:  off 

High-Frequency  Tesserals:  off 

Third -Body : 

off 

J22  :  off 

Central  Body 

J2 /M -Daily:  off 

These  models  were  chosen  to  simulate  the  SDP4  dynamics  described  in  [23] 


and  [24].  With  this  dynamic  modelling,  the  "SDP4-simulated"  DC  of  20  days 
of  observations  for  NSSC  13964  resulted  in  the  initial  conditions  shown  in 
Table  22.  Note  the  final  values  for  mean  semimajor  axis  in  Table  22  and 
Table  16  differ  by  about  75  meters.  Note  also  the  high  value  of  the 
weighted  RMS. 


Table 

22.  Results  of 

the  “SDP4-Simulated"  DC  for 

NSSC  13964 

Epoch 

Jan  1,  1985,  23  hours,  38  min 

17.903  sec 

Element 

A  priori  Value 

Final  Value 

Standard 

Deviation 

a  (km ) 

26554.3000 

26554.3768 

9 . 74558E-04 

e 

0.71 790000 

0.71 7871 78 

1  .  76389E-06 

i(°) 

63.0900000 

63.0848103 

9.71 891 E-02 

fi  (°> 

163.313000 

163.312091 

7.27630E-02 

M  (•> 

279.650000 

279.652677 

4.48756E-02 

M(°) 

12.6270000 

12.6445626 

1  .67609E-02 

Weighted 

RMS:  3.9805 

Mean  of  1950.0  Coordinates 

The  final  values  shown  in  Table  22  were  then  used  to  create  a  90-day 
prediction  file  using  the  dynamic  modelling  of  Table  20.  The  first  20 
days  of  prediction  (corresponding  to  the  DC  fit  interval)  were  compared  to 
the  fully-modelled  SST  predictions.  The  position  and  velocity  RMS  values 
of  the  errors  in  the  radial,  cross-track,  and  along-track  directions  are 


given  in  Table  23 
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Table  23.  Comparison  of  "SDP4-Simulated"  Predictions  and 
Fully-Modelled  SST  Predictions  for  NSSC  13964 
from  850102  to  850122 


Direction 

Radial 
Cross -Track 
Along-Track 
Total 


Position  RMS  (km)  Velocity  RMS  (km/sec) 


0.151 96D+01 
0.401 88D+01 
0 . 25957D+01 
0. 501 98D+01 


0.33271D-03 
0.401 20D-03 
0.171 27D-03 
0. 54862D-03 


An  effort  was  then  made  to  use  the  fully-modelled  version  of  SST  with  the 
shorter  fit  interval  (20  days)  to  generate  initial  conditions  for  the  SST 
predictions.  These  predictions  were  then  compared  to  the  predictions  that 
resulted  from  the  initialization  using  the  90-day  fit  interval.  The 
modelling  of  Tables  14  and  15  were  used  along  with  the  converged  value  for 
the  solar  radiation  coefficient  Cr  from  Table  16.  The  same  a  priori 
values  of  the  Keplerian  elements  as  Table  16  were  used.  The  results  are 
given  in  Table  24. 


Table  24.  Results  of  SST  DC  with  20-Day  Fit  for 
NSSC  13964 

Epoch:  Jan  1,  1985,  23  hours,  38  min,  17.903  sec 


Element  A  priori  Value  Final  Value 


a  (km) 


26554.3000 
0.71 790000 
63.0900000 
163.313000 
279.650000 
12.6270000 


Weighted  RMS:  1  .01 23 


Mean  of  1950.0  Coordinates 


26554.3030 

0.71791564 

63.0898234 

163.314403 

279.653454 

12.6279070 


Standard 

Deviation 

2 . 24820E-04 
4 . 24458E-07 
2.4171 9E-02 
1 . 79993E-02 
1  .1 1623E-02 
4 . 05857E-02 


The  weighted  RMS  of  1.0123  from  Table  24  compares  very  favorably  with  the 


weighted  RMS  of  0.92677  from  Table  16.  Note  also  that  the  final  values 
for  semimajor  axis  from  Tables  24  and  16  differ  by  less  than  1  meter. 
These  results  indicate  that  satisfactory  initial  conditions  can  be 
determined  with  the  shorter  interval  of  observations. 

The  final  values  of  the  Keplerian  elements  from  Table  24  and  the  same 
value  of  Cr  were  used  to  generate  a  prediction  file.  The  first  20  days 
of  predictions  (corresponding  to  the  DC  fit  interval)  were  compared  to  the 
predictions  that  resulted  from  the  fully-modelled  DC  with  the  90-day  fit 
interval  (Table  16).  The  position  and  velocity  RMS  values  of  the  errors 
in  the  radial,  cross-track,  and  along-track  directions  are  given  in  Table 
25. 


Table  25.  Comparison  of  SST  Predictions  (Using  20-Day  Fit 
Interval)  and  SST  Predictions  (Using  90-Day  Fit 
Interval)  for  NSSC  13964  from  850102  to  850122 

Direction 

Position  RMS  (km) 

Velocity  RMS  (km/sec) 

Radial 

Cross -Track 

Along -Track 

Total 

0.45229D-01 

0.1 0683D+01 

0.1 4532D+00 
0.10791D+01 

0 . 1 0455D-04 

0.1 1 162D-03 

0.11 533D-04 

0.1 1270D-03 

A  comparison  of  Table  23  and  Table  25  demonstrates  the  effects  of  reduced 
dynamic  modelling  for  highly  eccentric  orbits.  The  enhanced  dynamic 
modelling  in  SST  resulted  in  much  smaller  errors. 
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The  next  effort  was  to  evaluate  the  predictions  beyond  the  fit  interval. 
When  the  final  values  of  Table  24  were  used  to  generate  a  prediction  in 
the  850122  to  850212  interval,  results  were  very  good,  as  seen  in  Table 
26. 


Table  26.  Comparison  of  SST  Predictions  (Using  20-Day  Fit 
Interval)  and  SST  Predictions  (Using  90-Day  Fit 
Interval)  for  NSSC  13964  from  850122  to  850212 

Direction 

Position  RMS  (km) 

Velocity  RMS  (km/sec) 

Radial 

Cross -Track 
Along-Track 

Total 

0.1 3058D+00 

0.11 305D+01 

0. 73424D-01 

0.11 404D+01 

0.1 7235D-04 

0 . 78288D-04 

0.59191 D-05 

0 . 80380D-04 

A  comparison  of  the  radial  and  a  long-track  errors  in  Tables  25  and  26 
demonstrates  that  position  errors  beyond  the  fit  interval  very  closely 
match  the  position  errors  calculated  for  the  fit  interval. 

Figures  59  through  61  show  the  radial,  cross-track,  and  along-track  errors 
for  the  predictions  during  the  DC  fit  interval.  Figures  62  through  64 
show  the  prediction  errors  for  the  following  20  days.  Each  of  these 
figures  show  two  apparent  curves  because  the  1/4  day  comparison  interval 
is  approximately  commensurate  with  the  1 2-hr  orbital  period. 
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Figure  63.  NSSC  13964  Cross-Track  Differences 
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CHAPTER  4 


CONCLUSIONS  AND  RECOMMENDATIONS 


4.1  CONCLUSIONS 

The  primary  objective  of  this  thesis  was  to  test  the  semianalytical  satel¬ 
lite  theoryf;as  implemented  in  the  Charles  Stark  Draper  Laboratory  version 
of  the  Goddard  Trajectory  Determination  System  against  long  arcs  of  real 
data  of  highly  eccentric  orbits.  The  real  data  consisted  of  a  9-year  arc 
of  NORAD  element  sets  for  NSSC  9829,  the  Soviet  Molniya  2-17  spacecraft,* 
an  18-month  arc  of  element  sets  for  NSSC  14095,  the  ESA  Exosat/  and  an 
18-month  arc  of  NORAD  observations  and  element  sets  for  NSSC  13964,  the 
Soviet  Molniya  1-57  spacecraft.  .. 

An  introductory  chapter  reviewed  the  uses  of  elliptical  orbits  and  the 

i  *  •  i 

need  for  long-term  predictions^  It  also  reviewed  'the  application  of 
artificial  satellite  theory  to  elliptical  orbits,  including  a  summary  of 
the  previous  applications  of  satellite  theory  to  highly  elliptical  orbits. 

The  overall  test  methodology  involved  four  main  steps: 


1 .  Pre-processing  NORAD  data 

2,  Determining  an  initial  state  using  an  SST  DC 


"•  wVvVVVvV' 
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3i  Generating  an  SST  prediction  file 
4 .  Evaluating  the  SST  prediction 

For  NSSC  9829,  the  Soviet  Molniya  2-17,  the  usage  of  element  sets  resulted 
in  an  initial  vehicle  state  with  a  weighted  RMS  of  0.2893  (Table  4).  Best 
results  were  obtained  using  a  6x6  field  with  GEM  9  coefficients.  This 
initial  state  was  propagated  about  2330  days  in  an  SST  process  to  create  a 
9-year  arc  of  predictions  which  were  compared  to  the  NORAD  mean  elements. 
The  slowly  varying  Keplerian  element  histories  and  the  mean  anomaly 
difference  plot  clearly  demonstrate  the  ability  of  the  semianalytic  theory 
to  predict  the  dominant  motion.  Difference  plots  for  semimajor  axis 
exhibited  a  frequency  very  similar  to  the  semimajor  axis  resonance  period. 
Efforts  to  include  an  8x8  field  with  GEM  9  coefficients  resulted  in  a  very 
unsatisfactory  fit  (a  weighted  RMS  of  2.071).  Exhaustive  testing  then 
validated  the  SST  models  with  the  8x8  field. 

Actual  observation  data  was  available  for  NSSC  13964,  the  Soviet  Molniya 
1-57.  This  data  was  pre-processed  to  determine  the  initial  state  (includ¬ 
ing  the  solar  radiation  coefficient  Cr )  and  resulted  in  a  much  improved 
weighted  RMS  of  0.92677  (Table  16,  with  the  uncertainties  of  Tables  12  and 
13).  These  initial  conditions  were  used  to  generate  an  SST  prediction  of 
about  250  days.  Comparison  and  difference  plots  exhibited  a  significant 
improvement  over  the  NSSC  9829  case  as  shown  in  Table  27. 


Table  27 

Comparison  of  Prediction  Results 

Using  Element  Sets 

and  Observations 

Element 

NSSC  9829 

NSSC  13964 

Mean  Difference 

Mean  Difference 

a(m) 

250.46 

67.695 

e 

0.22437E-03 

0. 76003E-04 

i(°) 

0 . 33941 E-01 

0. 78965E-02 

ft  (°) 

0.22964 

0.1 5868E-01 

0)  (°) 

0.11699 

0.1 7679E-01 

M<°) 

1 .0960 

0.19882 

The  force  models  used  in  SDP4,  actual  observations,  and  short  fit  span 

were  included  in  an  SST  DC  to  create  "SDP4-simulated"  predictions.  It  was 
demonstrated  that  much  shorter  intervals  of  observations  (20  days)  could 
be  used  to  initialize  the  SST  predictions. 

For  the  Exosat  orbit  the  mean  difference  between  the  NORAD  and  SST 

predictions  for  semimajor  axis  was  approximately  77  kilometers  over  the 

250  days  of  the  comparison  interval.  The  mean  differences  for  the 
eccentricity,  inclination,  and  argument  of  perigee  were  0.0013,  0.102°, 

and  0.124°,  respectively. 

This  study  -has  shown  that  the  semianalytical  theory  results  in  very 

satisfactory  predictions  for  orbits  of  high  eccentricity.  For  Molniya- 
type  orbits,  the  results  obtained  when  using  actual  observational  data 
were  superior  to  the  results  obtained  when  using  element  sets.  For  orbits 
of  even  higher  eccentricity,  the  SST  results  were  still  very  satisfactory. 
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4.2  FUTURE  WORK 

The  use  of  actual  NORAD  observations  significantly  improved  the  weighted 
RMS  of  the  differential  correction.  But  there  were  only  about  160  days  of 
element  data  beyond  the  fit  interval  for  NSSC  13964  to  make  comparisons. 
A  longer  arc  (through  1986)  of  element  sets  would  permit  more  meaningful 
comparisons.  Specifically,  absence  of  secular  drift  in  a  difference  plot 
for  the  semimajor  axis  would  confirm  that  the  8x8  tesseral  field  was 
sufficient  for  long-term  predictions  of  Molniya  orbits.  Difference  plots 
for  a  longer  interval  would  also  determine  if  the  trends  in  longitude  of 
ascending  node  and  argument  of  perigee  (Figures  54  and  55)  were  secular  or 
periodic. 

The  orbits  of  the  spacecraft  in  this  study  were  chosen  so  that  the  perigee 
remained  outside  the  main  portion  of  the  atmosphere.  Thus  the  drag 
effects  were  very  slight  and  the  period  remained  constant.  This  enabled  a 
single  tesseral  resonance  modelling  to  remain  valid  over  the  entire 
prediction  span. 

For  some  orbits,  the  lunar-solar  perturbations  hold  the  perigee  height  in 
the  200  to  350  km  range.  (An  example  is  NSSC  10723,  the  International 
Ultraviolet  Explorer  rocket  body.)  The  orbit  is  being  continually  per¬ 
turbed  by  drag  and  the  orbital  period  decreases  significantly  over  long 
periods  of  time.  Thus  the  software  modelling  for  the  tesseral  geopoten¬ 
tial  harmonic  resonance  would  have  to  be  variable  over  the  prediction 
interval,  or  several  shorter  intervals  could  be  carefully  "piece-mealed" 
together. 
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The  accurate  calculation  of  drag  effects  for  this  type  of  orbit  is  also 


complicated  by  the  long-term  changes  in  atmospheric  density.  Software 
modelling  of  drag  would  probably  require  a  "smoothed"  modelling  of  density 
over  the  solar  cycle.  Long  period  variations  in  both  the  solar  activity 
and  the  geomagnetic  index  should  be  considered. 


Extensions  to  SDP4 

analysis 

to  tie 

down  the 

particular 

source 

of 

mismodelling  which 

makes  the 

major 

contribution 

to  error 

should 

be 

considered . 

The  analysis  of  the  Exosat  orbit  was  based  on  real  data  in  the  form  of 
NORAD  element  sets.  Preliminary  analysis  of  Exosat  using  observation  data 
has  proven  unsatisfactory  probably  due  to  the  sparseness  of  the  data. 
Long  arcs  of  ESA  observation  data  (with  the  appropriate  error  statistics 
and  station  locations)  could  be  combined  with  NORAD  observational  data. 
This  would  probably  provide  more  accuracy  in  the  same  manner  that  NORAD 
observations  improved  the  Molniya  predictions. 

Alternatively,  we  note  that  our  interest  in  the  EXOSAT  was  motivated  by 
the  desire  to  test  the  SST  for  an  orbit  with  an  eccentricity  much  greater 
than  a  Molniya.  However,  the  extreme  altitude  associated  with  most  of  the 
EXOSAT  orbit  limits  the  available  observation  data.  Perhaps  it  would  be 
desirable  to  choose  an  orbit  with  a  smaller  maximum  altitude  for  this 
test.  The  orbits  of  Elektron  2  and  4  (e  =»  .8),  COS-B  (e  K  .85),  and 
ISEE1/2  (e  “  .91)  should  be  considered  in  this  context. 
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APPENDIX  A 


SOURCE  CODE  FOR  UTILITIES 


The  following  pages  contain  the  IBM  FORTRAN  source  code  for  the  utilities 
used  in  this  study.  Listings  are  available  for: 


ADCEDIT 

RUNADC 

SDP4 

RDORBl 

PLOTTER 

RUNADCOB 

ASTRON 

PRENUT 
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PROGRAM 


ADCEDIT 


******** 

FUNCTION 

******** 


This  program  creates  a  control  deck  and  a  data  deck  to  be  used  in 
a  CSDL  PL0T4B  plotting  program.  It  can  be  used  to  visually  edit 
NORAD  bevel  vector  data. 


INPUT . . . . 

FILE  1  NORAD  bevel  vector  element  sets. 

OUTPUT  »*■•*******••****••**•*•*•*•**•••.*•*.••*•* 

FILE  40  PL0T4B  control  deck. 

FILE  41  PLOT 4B  data  set. 

SUBROUTINES  CALLED  *••***»*•**•**********•»•*•**•• 

none 


USAGE 


File  1  data,  the  NORAD  bevel  vector  element  sets,  can  be  no  more 
than  2000  lines,  i.e.,  1000  element  sets.  After  creating  the  PLOT 
4B  control  deck,  change  all  occurrences  of  'YR'  to  the  year  of 
the  desired  graph. 

Invoke  the  PLOT  4B  program  by  entering  the  command 
'ADC . .  EDIT  .  PLOT  .  FORM  IDENT(  AOC . .  EDIT  .  PLOT  ) ' 

. * . HISTORY 


VERSION:  September  1986 

Fortran  program  for  the  IBM  3081  and  3033. 

ANALYSIS 
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C  Martin  E.  Fieger,  Capt,  USAF  AFIT  /  MIT 

C 

C  PROGRAMMER 

C  Martin  E.  Fieger,  Capt,  USAF  AFIT  /  MIT 

C 

C 

C 

C 

C 

C«»». .............  DECLARATIONS  . . . 

c 

c 

c 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C 

C 

C  Character  Variables  »»•»*•*»*•»*•••••«•».•«•••••*•» 

C 

C  none 

C 

C  Dimensions  •*•«*••»•»»*••••••»••»»«•»•»»•«»•»•*•••* 

C 

DIMENSION  ELEM  (10000,10) 

C 

C  Data  statements  •*»«•»••••*••*»•*•••»•.»•••»•«•.••• 


C 

C  none 
C 

C. .......... ......  begin  PROGRAM  . . 

C 

C  Set  counter  to  zero. 

C 

NUMOBS  ■  0 

C 

C  Read  In  a  new  mean  element  set. 

C 


100  CONTINUE 
C 

READ  ( 1 , 1000, END '200 )  NUMSV,  DATE,  XNDT2  ,  XNDD6, 

•  I  EXP ,  BSTAR ,  IBEXP,  XINCL,  XNODE ,  EO,  OMEGA, 

*  XMO ,  XNO 
C 

C  Calculate  real  values. 

C 

XNDD6  .  XNDD6  *  (10. DO  **  IEXP) 

BSTAR  .  BSTAR  *  (10. DO  **  IBEXP) 

C 

C  Increment  counter. 

C 

NUMOBS  *  NUMOBS  +  1 
C 

C  Store  data  in  array  ELEM. 

C 

ELEM  (NUMOBS. 1)  *  DATE 


( NUMOBS ,2 ) 
(NUM0SS.3) 
(NUMOBS, 4) 
(NUMOBS, 5) 
(NUMOBS, 6) 
(NUMOBS. 7) 
(NUMOBS, S) 
(NUMOBS, 9) 
(NUMOBS,  tO ) 


GO  TO  100 


PLOT 4B  OAT A  DECK 


Write  the  PL0T4B  date  deck. 


200  CONTINUE 


DO  300  I -=1, NUMOBS 


WRITE  (41,2000)  ELEM  (1,1),  ELEM  (1,2),  ELEM  (1,3),  ELEM  (1,4) 


300  CONTINUE 


DO  400  1*1, NUMOBS 


WRITE  (41,2000)  ELEM  (1.1),  ELEM  (1,5),  ELEM  (1,6),  ELEM  (1,7) 


400  CONTINUE 


DO  500  1*1. NUMOBS 


WRITE  (41,2000)  ELEM  (1,1),  ELEM  (1,8),  ELEM  (1,9),  ELEM  (1,10) 


500  CONTINUE 


PL0T4B  CONTROL  DECK 


Write  the  PL0T4B  control  cards. 


WRITE  (40,3000)  NUMOBS,  NUMSV 


WRITE  (40,4000) 


WRITE  (40.3000)  NUMOBS,  NUMSV 


WRITE  (40.5000) 


WRITE  (40,3000)  NUMOBS,  NUMSV 


WRITE  (40.6000) 


600  CONTINUE 

Write  final  P10T4B  control  card. 

WRITE  (40,7000) 

STOP 

*****************  FORMAT  STATEMENTS  . . . 

Input  cards  (  file  1  ). 

1000  FORMAT  (2X,I5,11X1F6.0,9X,D108,2(1X,D6.5,I2)/ 

*  7X  ,2( 1X.D8 .4) .  1X.D7  7,2( 1X.D8 .4)  .1X.D1 1 .8) 

Formatted  data  card 

2000  FORMAT  (  G12 . 5 ,3(2X ,614 . 7  )  ) 

PL0T4B  control  card 


3000  FORMAT 

•  (  ' “DAT  A  NUMVAR  4 

•  '  NOREWIND  '.14 

•  FMTDATA  1 

•  '(  G12. 5 ,3( 2X ,G14 . 7  )  ) 

•  '  ‘TITLE  MAINTITL  0 

•  ' SV ' , 15 ■ 

•  '  SUBTITL  0 

•  'NORAD  MEAN  ELEMENTS 

•  '  AXTITL  1 

•  'DATE  (YVDDD) 

•  ' 19YR 

•  ' *PLOTMOD  REMARK 


•  '  BADPOINT  1002 

•  ' *PLOT  TITLPLOT  0 

•  '  LEGEND  0 

•  '  2 

4000  FORMAT 

•  ( ' *PLOTMOD  NOZERO 

•  '  BLKPLOT 

•  ' XNDT2 
*'312 

•  '  YR001. 

•  '*PLOTMOD  NOZERO 

•  '  BLKPLOT 

•  ' XN0D6 

•  '  3  1  3 

•  '  YR001 . 

•  '  ‘PLOTMOD  NOZERO 


BLKPLOT 


5000  FORMAT 


( ' ‘PLOTMOD  NOZERO 


BLKPLOT 


'•PLOTMOD  NOZERO 


BLKPLOT 


'•PLOTMOD  NOZERO 


FORMAT 

('•PLOTMOD  NOZERO 
BLKPLOT 


•PLOTMOD  NOZERO 


'•PLOTMOD  NOZERO 


BLKPLOT 


7000  FORMAT  (  ' »END  LAST' 


■»» 


Vv  ’v  v  V  " 

> -V.VA  .VA 


mm  mm  *.  k  »v».  r.-*  j 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  RUNADC 


FUNCTION 


This  program  pre-processes  NORAO/ADCOM  general -perturoat ion 
mean  orbital  elements  for  use  in  a  GTOS  run. 


INPUT  . . . . . . . 

FUE  1  NORAC  bevel  vector  element  sets,  e  g  , 

ADC . .  ELS  DATA 

where  .  is  the  satellite  designator 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OUTPUT  . . . . . . . 

FILE  2  GTOS  observation  cares,  e.g., 

ADC . . MEAN. OBS. DATA  or 

ADC . .OSCU.  OBS  .DATA 

FILE  3  NORAD  s i ngl e-averagee  elements,  e.g., 

ADC . .NORAD. MEAN. DATA  or 

ADC . .NORAD. OSCU. DATA 

FILE  €  Printer  messages. 

SUBROUTINES  CALLED  ••*»••**•••*•*•••*«*»»»»*»»•**»*»**»**»•»*»**»• 
SDP4  ADDTIM  CALPAK  JULIAN  EQUIN  KEPEQN 


HISTORY 


VERSION:  September  1986 

Fortran  program  for  the  IBM  3081  and  3033. 

ANALYSIS 

(unknown)  --  NORAD  /  ADCOM 

PROGRAMMER 
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( unknown ) 


NORAD  /  ADCOM 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c> 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


MOD  I F I E  C  .  Paul  d.  Cef  ola 

October  1981  -  Charles  Stark  Draper  Laboratory 


1.  Constructed  double-precision  version. 

2.  Modified  program  to  run  on  IBM  at  CSDL . 

MODIFIED  .  Martin  Fleger 

duly  1986  . AFIT  /  USAF 

1.  Pre-process  NORAD  mean  element  sets. 

2.  Write  NORAD  mean  element  sets  in  GTDS  observation  card 
format . 


MODIFIED  .  Lee  W.  Early,  dr. 

duly  1986  .  Charles  Stark  Draper  Laboratory 


1  .  Cleaned  up  code  . 

2.  Added  debugging  write  switch. 

MODIFIED  .  Martin  E  Fleger 

August  1986  .  AFIT  /  USAF 

1.  Modified  program  to  write  the  NORAD  elements  into 
a  temporary  data  set  to  be  used  for  plotting. 

2.  Added  switch  to  Subroutine  SDP4  to  output  either  mean 
data  (IUPSP*MEAN)  or  osculating  data  ( IUPSP*OSCU ) . 

. .  declarations  . . . 


IMPLICIT  double 

PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION 

EONELMt  6 )  ,KEPELM<6) 

Character  Variables  ***’ 

CHARACTER  *4  ISET, 

ITYPE, 

IUPSP 

Dimensions  ************' 

DIMENSION  PO S ( 3  ) 

,VEL( 3 ) 

/Cl/ 


COMMON  /Cl / 

CK2 

,CK4 

.  E6A 

.00MS2T 

S 

, TOTHRD 

,XJ3 

,  XKE 

XKMPER 

, XMNPDA 

.AE 

c 

c 


/  C2  ■ 


COMMON  / C2 / 


COMMON  /El/ 

XMO 

,  XNODEO 

.OMEGAO 

.EO 

* 

XINCL 

,  XNO 

, XNDT20 

, XNDD60 

* 

BSTAR 

.X 

,  Y 

,Z 

* 

XDOT 

DS50 

,  YDOT 

,  ZDOT 

.EPOCH 

DATA  STATEMENTS 


DATA  00 
DATA  SO 
DATA  X02 
DATA  X04 
DATA  GM 


120. 

78 

1 .082616 
1 .65597 
3.986008 


DATA  TIMTOL 


BEGIN  PROGRAM 


Read  control  canos 


Read  satellite  designator,  NORAD  generator  type, 
card  format,  type  of  short-period  motion  desired, 
start  date  of  data  files,  end  date  for  the 
observation  card  file,  end  date  for  the  element 
data  file,  search  Interval,  and  the  suppression 
Interval . 


READ  (5,1000)  INTLSV,  ISET,  ITYPE,  IUPSP,  DST  ART ,  ENDOBS, 
•  DSTOP ,  INTVL ■  IDE  LTA 


Error  on  control  card? 


IF  (ITYPE  .NE.  ’ TRNS ' )  THEN 
WRITE  (6,2000)  ITYPE 


STOP 
ENC  IF 


IF  (ISET  . NE  SDP4  )  THEN 
WRITE  (6,2010)  ISET 


’ 


'.I 


i  Ll  l.l  4.1 


Output  control  caras  to  printer 


WRITE  (6,1010)  INTLSV,  ISET,  ITYPE,  IUPSP,  DST ART ,  ENDOBS, 
•  DSTOP  ,  INTVL,  1DELTA 


INITIALIZE 


IF  (  ISET  .EQ.  'S0P4-  )  I EPT  =  3 


SET  FILE  PARAMETERS 


Compute  constants 


.500  *  XJ2  *  AE  * "2 


.37500  *  XJ4  *  AE’*4 


( (00  -  SOI  *  AE  /  XKMPER )  • >  4 
AE  *  (  1.00  +  SO  /  XKMPER  ) 


Convert  start  date,  end  obs  date,  stop  date 
to  internal  units  and  the  interval  periods 


to  seconds. 


DST  ART  •  1.06 


ENDOBS  *  1.06 


DSTOP  ■  1.06 


INTVL  '  86400.00 


IDELTA  *  3600. DO 


Write  first  GTDS  observation  ca”d  on  file 


WRITE  (2,4000)  'OBSCARD  ' 


Print  first  GTDS  observation  card. 


WRITE  (6,5000)  'OBSCARD 


SET  CONTROL  VARIABLES 


Set  computed  GO  TO  switch. 


A .Sv>v.vv.v-  ;>»v '  v  .  's  v  v 


s'.'.  V’  , 
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C  Set  first  element  set  switch. 

C 

NUMOES  =  0 

C 
C 
C 
C 

r  ..»*.*.*********».»** 

C  READ  MEAN  ELEMENT  SET 

r  »„•***********.**»»*. 

C 

C 

C 


c  ..........  READ  MEAN  ELEMENT  SET  .......... 

C 

C  Read  in  a  new  mean  element  set. 

C 

100  CONTINUE 

READ  ( 1 ,300C,END=500)  NUMSV,  IYRN,  IDOYN.  DAYPTN.  XNDT2N , 


*  XNDD6N ,  IEXPN,  BSTARN ,  IBExPN,  XINCLN,  XNDDEN ,  EDN, 

*  OMEGAN,  XMON,  XNON 

IF  (XNON  .LE.  0.00)  STOP 
C 

C  Check  satellite  designator. 

C 

IF  ( INTLSV  NE  NUMSV  )  THEN 
WRITE  (6,2020) 

STOP 

END  IF 

C 

C  Fine  epoch  of  new  mean  element  set. 

C 

E POCHN  *  IYRN  •  1000  +  IDOYN  *  DAYPTN 

SECS  *  86400. DO  *  DAYPTN 

YEAR  =  IYRN 

C 

C  Convert  format  of  epoch  to  YYMMDDHHMMSS . SSS 

C 

CALL  JULIAN  ( DAY JUL , SECJUL ,  YEAR , 1 . D0,0 . DO, 

*  0 . DO ,0 . DO ,0  DO ) 

DAYJUL  =  DAY JUL  +  IDOYN 

SECJUL  =  SECJUL  +  SECS 

CALL  CALPAK  (YMD.HMS,  DAY JUL , SECJUL ) 

EPOKEN  =  YMD  ■  1  ,D6  +  HMS 


C 

C  Print  epoch  t imes  . 

C 

WRITE  (6,5010)  EPOCHN ,  EPOKEN 
C 
C 
C 

C  ..........  ACCEPT  OR  REJECTS  ••••• 

C 

C  Branch  to  appropriate  part  of  program. 
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New  epoch  before  beginning  of  run' 
150  IF  ( EPOKEN  .IT.  DSTART )  GO  TO  100 


FIND  NEXT  EPOCH 


C  Update  old  mean  element  set  with  new  mean  element 

C  set . 

C 


200 


C 

C 

C 


C 

C 

C 


C 

C 

C 


CONTINUE 


EPOKEO 

= 

EPOKEN 

EPOCH 

= 

EPOCHN 

XNDT20 

= 

XNDT2N 

XNCDSO 

= 

XNDD6N 

I  EXP 

= 

IEXPN 

IBEXP 

= 

IBEXPN 

bstar 

s 

bstarn 

XINCL 

= 

XINCLN 

XNOOEC 

= 

XNOOEN 

EO 

r 

EON 

OMEGAO 

= 

OMEGAN 

XMO 

s 

XMON 

XNO 

s 

XNON 

Save  epoch  time. 


DAYSAV 

*  DAY JUL 

SECSAV 

=  SECJUL 

XJUL 

-  DAvsav 

*  SECSAV 

/  86400 . DO 

YMDSAV 

=  YMD 

HMSSAV 

=  HMS 

Increment 

epoch  by 

specified  search  interval 

CALL 

ADDTIM  (YMD 

,HMS, 

DAYJUL, SECJUL,  XINTVL, 

XDATE 

'  YMD  *  1 

.06  + 

HMS 

Restore  epoch  time 

DAY JUL 

=  DAYSAV 

SECJUL 

=  SECSAV 

YMD 

=  YMDSAV 

HMS 

«  HMSSAV 

TIMTOL  ) 


C 

C 


Increment  epoch  by  specified  suppression  interval. 


CALL  ADDTIM  (YMD.HMS,  DAYJUL , SECJUL ,  DELTA,  TIMTOL) 
XDELTA  =  YMD  *  1.D6  +  HMS 


Read  next  mean  element  set. 

2 


IGO 

GO  TO  100 


FIND  POSITION  AND  VELOCITY 


SET  ORBIT  GENERATOR 
PARAMETERS 


New  epoch  within  hourly  suppression  zone  of  previous 
epoch" 


300  IF  (EPOKEN  . LT .  XDELTA)  GO  TO  200 


Compute  orbit  generator  parameters. 


XNDD60 

r 

XNDD60  *  <10.00 

IEXP) 

XNODEO 

x 

XNODEO  *  DE2RA 

OMEGAO 

s 

OMEGAO  •  DE2RA 

XMO 

X 

XMO  •  DE2RA 

XINCL 

X 

XINCL  '  DE2RA 

TEMP 

c 

TWOPI  /  ( XMNPDA 

•  XMNPDA 

XNO 

= 

XNO  *  TEMP  *  XMNPDA 

XNDT20 

X 

XN0T20  •  TEMP 

XND060 

c 

XN0D60  *  TEMP  / 

XMNPDA 

Compute  more  orbit  generator  parameters 
«  ( XKE  /  XNO  )  *•  TOTHRD 

*  1.5D0  *  CK2  •<3.DO*(DCOS(XINCL)**2.DO)  -  1 . DO  1  / 
DO  -  E 0 * E 0  }  *'  1 . 5D0 ) 

■  TEMP  /  (At  *  A 1 ) 

-  A 1  *  ( 1 . DO  -  DELI  *  ( . 5D0  *  TOTHRD  * 

1  •  ( 1 . DO  ♦  (  134. DO/81  D0)*DEL1 )) 1 
«  TEMP  /  (AO  •  AO) 

*  XNO  /  ( 1 .DO  ♦  DELO  ) 

*  BSTAR  •  (10. DO  •*  IBEXP)  /  AE 

Set  input  parameters  for  call  to  orbit  generator 


TSINCE  *  0.00 

IFLAG  *  1 


•j  v:.*  v; 


c 
c 
c 

c  **» . .  EVALUATE  ORBIT  GENERATOR  ****** 

C 

C  Which  analytical  orbit  theory? 

C 

C  Note  that  SDP4  is  the  only  subroutine  that  has 

C  been  modified  to  include  the  additional  parameter 

C  IUPSP. 

C 

GOTO  (351,352,353,354.355),  IEPT 
C 

C  . SGP . 

C 

351  CALL  SGP  ( I  FLAG  .TSINCE  ) 

GO  TO  360 

C 

C  . SGP4  . 

C 

352  CALL  SGP4  ( I  FLAG .TSINCE ) 

GO  TO  360 

C 

C  S0P4  . 

C 

353  CALL  SDP4  ( IFLAG, TSINCE , IUPSP ) 

GO  TO  360 

C 

C  SGP8  . 

C 

354  CALL  SGP8  ( I  FLAG , TSINCE  ) 

GO  TO  360 

C 

C  SDP8 . 

C 

355  CALL  SDP8  ( I F LAG  .TSINCE  1 

C 

C 

c  WRITE  OUTPUT 

C 

C  Compute  position  in  kilometers. 

360  X  *  X  *  XKMPER  /  AE 

V  *  Y  *  XKMPER  /  AE 

Z  *  Z  *  XKMPER  /  AE 

C 

C  Compute  velocity  In  kilometers/second. 

C 

C 

Ox  *  XOOT 

DY  *  YOOT 

DZ  *  ZDOT 

C 


DX 


DX  *  (XKMPER  /  AE)  *  (XMNPDA  /  86400. DO) 


DY 

02 


DY  *  (XKMPER  /  AE)  ■  (XMNPDA  /  86400.00) 
02  *  (XKMPER  /  AE)  *  (XMNPDA  /  86400. DO) 


Write  position  on  GTDS  oPservation  card  file. 

IF  (  EPOKEO  . LE.  ENDOBS  )  THEN 

WRITE  (2,4010)  21  .EPOKEO, X.X 

WRITE  (2,4010)  22,  EPOKEO, Y,Y 

WRITE  (2,4010)  23, EPOKEO, 2, 2 

Write  velocity  on  GTDS  observation  card  file. 

WRITE  (2,4010)  24, EPOKEO, OX, OX 

WRITE  (2,4010)  25,  EPOKEO, DY.DY 

WRITE  (2,4010)  26, EPOKEO, 02,02 


Print  position  in  GTDS  oPservation  card  format. 

WRITE  (6,5020)  2  1  ,  EPOKEO  ,  X , X 
WRITE  (6,5020)  22  ,  EPOKEO , Y , Y 
WRITE  (6,5020)  23  , EPOKEO . 2 , 2 

Print  velocity  in  GTDS  observation  card  format. 

WRITE  (6,5020)  24  ,  EPOKEO , Ox , OX 
WRITE  (6,5020)  25  . EPOKEO ,DY ,DY 
WRITE  (6,5020)  26  ,  EPOKEO  ,02 ,02 

Convert  to  equinoctial  elements. 

POS(  1  )  *  X 
POS( 2 )  =  Y 
POS  (  3 )  *  2 


CALL  EQUIN  (  EONE LM ,  RETRO ,  POS.VEL,  GM,  .TRUE.  ) 

SMA  «  EQNELM(I) 

XH  *  EQNELM(2) 

XK  x  EONE  LM( 3 ) 

XP  *  EQNELMU) 

XO  *  EQNELM( 5 ) 

XML  *  EONE  LM( 6 ) 


c 

c 


Convert  to  Keplerlan  elements. 


CALL  KEPEQN  (  KEPELM,  EONELM, RETRO  ) 

C 

ECC 
XINC 
XLAN 
AP 
XMA 
C 
C 
C 

PR 
APR 
C 
C 
C 

XML 
XINC 
XLAN 
AP 
XMA 
C 
C 
C 

WRITE  (3,6000)  XJUL  ,  EPOKEO , X , Y . Z ,0X ,DY ,DZ , SMA , PR , APR , ECC , 

•  X  INC, XLAN, AP, XMA, XH.XK.XP.XQ, XML 

C 

C  Increment  the  counter. 

C 

NUMOBS  *  NUMOBS  +  1 
C 

C  ........ 

C  END  LOOP 

C  . 

c 

c 

c 

C  New  epoch  beyond  end  of  run? 

C 

400  IF  (EPOKEN  . GT .  DSTOP  )  GO  TO  500 
C 

C  New  epoch  beyond  weekly  suppression  zone  of  previous 

C  epoch? 

C 

IF  (EPOKEN  . GT .  XDATE)  GO  TO  200 
C 


=  KEPELMt  2 ) 

=  KEPELM(  3 ) 

*  KEP£LM(4) 

=  K£PELM( 5 ) 

*  KEPELM(6) 

Calculate  perigee  radius  and  apogee  radius. 

=  SMA  •  (  1 .DO  -  ECC  ) 

=  SMA  »  (  1 .DO  +  ECC  ) 

Convert  angular  quantities  to  degrees. 

=  XML  /  DE2RA 

=  XINC  /  DE2RA 

=  XLAN  /  0E2RA 

=  AP  /  DE2RA 

=  XMA  /  DE2RA 

Write  days  elapsed  ana  orbital  elements. 


IGO  '  3 

GO  TO  100 


Read  next  mean  element  set 


END  RUN 


C 

1000  FORMAT  ( 33X , 1 5/3 ( 3 1 X . A4/ ) , 3( 3 1 X , F7 . 0/  )  ,3 1 X  ,  12/3 IX  ,  12  ) 

C 

1010  FORMAT  (IX, 'SV' , 1 5 ,2X ,3(  A4 ,2X ) ,3( F7 . 0 ,2X  )  ,2( 12 ,2X  )  ) 

C 

C  Error  messages. 

C 

2000  FORMAT (  '  A  card  type  of  ',A4,'  is  illegal.'  / 

*  '  You  must  use  a  transmission  card  to  input  data.') 

2010  FORMAT (  '  Ephemens  type  ',A4,'  not  legal;  will  skip  this  case.') 
C 

2020  FORMAT (  '  Satellite  designators  do  not  match.'  ) 

C 

C 

C 

C  Input  cards  (  file  5  ). 

C 

3000  FORMAT  ( 2X , 15 , 1 1 X , 12 . 1 3 ,09 . 8 , 1 X ,D 10 . 8  ,2<  IX  ,D6 . 5 , 1 2  ) / 

*  7X  ,2( IX  ,DS .4 ) , IX ,07 .7 ,2( 1X.D8 . 4 ) , IX ,D1 1 . 8 ) 

C 

C  GTDS  observation  cards  (  file  2  ). 

C 

4000  FORMAT  (  A8  ) 

4010  FORMAT (  8X ,  13, 6X,  G21.15,  2G21.14  ) 

C 

C  Debugging  print  (  file  6  ). 

C 

5000  FORMAT (  IX.  A8  ) 

5010  FORMATf  '  NORAD:  '.G19.13.5X,  'GTDS:  '.G21.15  ) 

5020  FORMAT (  9X  ,  13. 5X,  G21.15.5X,  G21.14.5X,  G21.14  ) 

C 

C  Temporary  data  set  (  file  3  ). 

C 

6000  FORMAT ( 'DAY' ,G19. 10.8X.G21 . 15  / 
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'X  ' ,G19. 10.7X, 'Y  ,G19 . 10, 7X,  1  '.G19.10/ 

OX  ,G19. 10.7X, 'DY  ,G19. 10.7X. 'DZ  '.G19.10/ 

' SMA  ,G19. 10.7X, 'PR  ,G19 . 10.7X, ' APR' ,G19 . 10  / 

'ECC  ,G19 . 10.7X , ' INC  ,G19 . 10 , 7X , ' LAN ' , G 19 . 10  / 

AP  ,G19. 10.7X, 'MA  • ,G19. 10.7X, 'H  '.G19.10/ 

'K  ' ,G19. 10.7X, 'P  '  ,G19 . 10.7X. 'Q  '.G19.10/ 

'ML  '.G19.10  ) 

END 

BLOCK  DATA 

PURPOSE 

INITIALIZE  PORTION  OF  COMMON  BLOCK  /Cl/  USED  IN  THE 
NORAD/ADCOM  GENERAL  PERTURBATION  SATELLITE  THEORY 
PACKAGE . 


DIMEN-  LOCA- 

SION  TION  DESCRIPTION 


E6A 

TOTHRD 

XKE 

XKMPER 
XMNPDA 
AE 
X  J3 


10* •  (  -6  ) 

2  /  3 

GRAVIATIONAL  CONSTANT , ( ER/MIN)*  *3/2 
KM  PER  EARTH  RADIUS 
TIME  UNITS  PER  DAY 

EQUATORIAL  RADIUS  OF  EARTH ,( D . U . /E . R . ) 
NUMERICAL  value  OF  J3  ZONAL  HARMONIC 


VERSION  OF  OCTOBER  1981 

FORTRAN  SUBROUTINE  FOR  THE  AMDAHL  470/V8  AND  THE  IBM  3033. 


ANALYSIS 

P.  CEFOLA 

PROGRAMMER 

P.  CEFOLA 


--  CHARLES  STARK  DRAPER  LABORATORY 


CHARLES  STARK  DRAPER  LABORATORY 


DECLARATIONS 


IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON/C  1 /CK2 .CK4.E6A .OOMS2T , S , TOTHRD , 
*  XJ3.XKE.XKMPER. XMNPDA, AE 


DATA 

E6A 

/ 

1.0-6 

/ 

DATA 

TOTHRD 

/ 

.6666666700 

/ 

DATA 

XKE 

/ 

. 743669 16  ID  -  1 

/ 

£ 


DATA 

XKMPER 

/ 

6378  135D0 

/ 

DATA 

XMNPDA 

/ 

1440.00 

/ 

DATA 

AE 

/ 

1  .DO 

/ 

DATA 

XJ3 

/ 

- .253881D-5 

/ 

END 

BLOCK 

DATA 

PURPOSE 

C  INITIALIZE  COMMON  BLOCK  /C2/.  THIS  BLOCK  DATA  IS  FOR  USE 

C  WITH  THE  NORAD/ADCOM  GENERAL  PERTURBATION  PACKAGE. 

c 

c 

C  VAR  I  -  OIMEN-  LOCA- 

C  ABLE  SION  TION  DESCRIPTION 

C  . 

c 

C  descriptive  heading 

C 

C  DE2RA  DEGREES  TO  RADIANS  CONVERSION 

C  PI  MATHEMATICAL  CONSTANT  PI 

C  PI02  PI  /  2 

C  TWOPI  2  *  PI 

C  X3PI02  3  *  PI  /  2 

C 

C  . 

C 

C  VERSION  OF  OCTOBER  1981 

C  FORTRAN  SUBROUTINE  FOR  THE  AMDAHL  47Q/V8  AND  THE  IBM  3033. 

C 

C  ANALYSIS 

C  P.  CEFOLA  --  CHARLES  STARK  DRAPER  LABORATORY 

C 

C  PROGRAMMER 

C  P.  CEFOLA  --  CHARLES  STARK  DRAPER  LABORATORY 

C 
C 
C 

. . *  DECLARATIONS  . . * . 

C 

IMPLICIT  REAL  *8  (A-H.O -2) 

C 

C0MM0N/C2/DE2RA ,PI , P 102 .TWOPI . X3P 102 


C 

DATA  DE2RA  /  .1745329250-1  / 
DATA  PI  /  3. 141 59265DO  / 
DATA  PI02  /  1 . 5707963300  / 
DATA  TWOPI  /  €.2831853000  / 
DATA  X3PI02  /  4 . 7 1238898DO  / 


SUBROUTINE  $DP4(IFLAG,TSINCE .IUPSP) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PURPOSE 

SDP4  IS  THE  12  HR,  HIGH  ECCENTRICITY  AND  24  HR 
GEOSYNCHRONOUS  NORAO/ADCOM  GENERAL  PERURBATION  THEORY. 
THIS  SUBROUTINE  IS  TAKEN  FROM  THE  PROJECT  SPACETRACK 
REPORT  NO. 3,  DECEMBER  1980,  'MODELS  FOR  THE 
PROPAGATION  OF  NORAD  ELEMENT  SETS",  F.R. HOOTS  AND 
R.  L.  ROEHRICH. 


METHOD 

BROUWER  -  HORI  -  LANE 


CALLING  SEQUENCE 

SDP4( IFLAG.TSINCE ) 


PHYSICAL  PARAMETERS  USED  IN  THE  COMMENTS 


parml 

parm2 

parm3 

parm4 


meaning 
meaning 
meaning 
mean i ng 


SUBROUTINES  CALLED 

DEEP  FM0D2P  ACTAN 


ADC0M/D06  VERSION  OF  DECEMBER  1980 

FORTRAN  SUBROUTINE  FOR  THE  AMDAHL  470/V8  AND  THE  IBM  3033. 
THIS  IS  A  DOUBLE  PRECISION  VERSION  CONSTRUCTED  BY  P.  CEFOLA 


ANALYSIS 

ANALYST  --  P.  CEFOLA 

PROGRAMMER 

PROGRAMMER  --  P.  CEFOLA 


MODIFICATIONS 

Version  of  September  1986 

Martin  E.  Fleger  -  Captain,  USAF/AFIT/MIT 

Added  a  subroutine  parameter  IUPSP  to  enable  selection  of 
a  mean  output  without  the  update  for  short  perlodlcs 
(IUPSP'1),  or  an  osculating  output  with  the  update  for 
short  perlodlcs  (IUPSP'O). 
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c, ...... ..........  DECLARATIONS 

c 


c 

c 

c 


c 


c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


IMPLICIT  REAL*8  (A-H.O-Z) 

CHARACTER*4  IUPSP 

COMMON/E 1/XMO , XNODEO .OMEGAO , EO , XINCL , XNO , XNDT20  , 

1  XNDD60 .BSTAR ,X ,Y , Z ,XDOT , YOOT .ZDOT . EPOCH .DS50 

COMMON/C 1 /CK2 ,CK4 , E6A .00MS2T . S . TOTHRD , 

1  X03,XKE,XKMPER,XMNPDA.AE 

.............  BEGIN  PROGRAM  . . 

IF  ( I  FLAG  .EO.  0)  GO  TO  100 

Recover  original  mean  motion  (XNODP)  and  semimajor 
axis  (AODP)  from  input  elements. 

A 1  *  (XKE/XNO )  **TOTHRD 

COSIO  *  DCOS( XINCL) 

THETA2  *  COSIO  *  COSIO 

X3THM1  *  3.0D0  *  THETA2  -  1 .ODO 

EOSO  '  EO  *  EO 

BETA02  *  1 . DO  -  EOSQ 

BETAO  *  DSQRT ( BET A02 ) 

DELI  *  1.5D0  *  CK2  »  X3THM1  /  (A1  *  A1  *  BETAO  *  BETA02  ) 

AO  *  A1  *  (  1. DO-DELI  *  ( . 5D0  *  TOTHRD  +  DE L 1  * ( 1 . D0+ 134 . DO/8 1 . DO* 
1  DELI  ))) 

DELO  *  1.5D0  *  CK2  *  X3THM1  /  (AO  *  AO  *  BETAO  *  BETA02) 

XNODP  *  XNO  /  (  1 . DO  +  DELO  ) 

AODP  *  AO  /  (  1 . DO  -  DELO  ) 

Initial Izatlon 


For  perigee  below  156  km,  the  values  of 
S  and  00MS2T  are  altered. 


Q0MS24  *  Q0MS2T 


PERIGE  *  (  AODP  *  (  1.DO-EO)  -  AE)  •  XKMPER 

IF  (PERIGE. GE. 158. DO)  GO  TO  10 

S4  *  PERIGE  -  78.00 

IF  (  PERIGE. GT.  98. DO  )  GO  TO  9 

S4  *  20. DO 

9  Q0MS24  =  ((120. DO  -  S4)  *  AE  /  XKMPER  )• *4 . DO 
S4  «  S4  /  XKMPER  +  AE 

10  PINVSO  *  1 . DO  /  (AODP  *  AODP  *  BETA02  *  BETA02) 

SING  *  DSIN(OMEGAO) 

COSG  =  DCOS(OMEGAO) 

TSI  *  1 .DO  /  (AODP  -  S4) 

ETA  =  AODP  •  EO  *  TSI 
ETASQ  -  ETA  •  ETA 

EETA  *  EO  «  ETA 

PSISQ  «  DABS  ( 1 .DO  -  ETASQ) 

COE F  *  Q0MS24  *  TSI**4,D0 
COE  FI  «  COEF  /  PSISQ**3.5DO 

C2  *  C0EF1*XN0DP*(A0DP*( 1 .DO  +  1.5D0  *  ETASQ  +  EETA  * 

*  (4.00  +  ETASQ))  +  . 75DO  *  CK2  *  TSI  /  PSISQ  * 

*  X3THM1  *  (  8. DO  +  3. DO  *  ETASQ  *  (  8. DO  +  ETASQ))) 

Cl  *  BSTAR  *  C2 
SINIO  *  DSIN(XINCL) 

A30VK2  *  -  X03  /  CK2  *  AE  **3.D0 
X1MTH2  «  1 .DO  -  THETA2 

C4  *  2. DO  *  XNODP  *  COE  FI  *  AODP  *  BETA02  *  (  ETA* 

*  (2. DO  +  .500  *  ETASQ  )  ♦  EO  *  ( . 5D0  +  2. DO  *  ETASQ) 

*  -2. DO  *  CK2  *  TSI  / 

*  (AODP* PSISQ)* ( -3 .DO  *X3THM1*( 1 .DO  -2. DO  *  EETA  +  ETASQ 

*  ( 1 . 5D0- . 5D0*EETA ) )+ . 7500*X 1MTH2  *  (2. DO  *  ETASQ  -  EETA 


( 1 . DO  +  ETASQ))  *  DC0S(2.D0  *  OMEGAO))) 


C 

THETA4  *  THETA2  *  THETA2 
C 

TEMPI  =  3.00  *  CK2  *  PINVSQ  *  XNOOP 
C 

TEMP2  =  TEMPI  *  CK2  *  PINVSQ 
C 

TEMP3  *  1.25D0  *  CK4  *  PINVSQ  *  PINVSQ  •  XNODP 
C 

XMDOT  =  XNODP  +  .500  *  TEMPI  •  BETAO  *  X3THM1 

*  +  .062500  »  TEMP2  *  BETAO  * 

*  (  13.00  -  78.00  *  THETA2  +  137. DO  *  THETA4  ) 

C 

X1M5TH  *  1.00  -  5.00  •  THETA2 

C 

OMGDOT  =  -.500  *  TEMPI  *  X1M5TH  +  .062500  *  TEMP2  * 

*  (7. DO  -  114. DO  *  THETA2  + 

*  395. DO  *  THETA4  )  +  TEMP3  *  (  3 . DO  -  36 . DO  *  THETA2  + 

*  49. DO  *  THETA4  ) 

C 


! 

c 

c 

c 

c 

c 

c 

c 

c 


XHD0T1  *  -  TEMPI  •  COSIO 

X NODOT  «  XHD0T1  +  (.5D0  *  TEMP2  *  (  4.00  -  19.00  •  THETA2  ) 

*  +  2.00  •  TEMP3  •  (  3.00  -  7.00  •  THETA2  ))  *  COSIO 

XNODCF  *  3.500  *  BETA02  *  XHD0T1  *  Cl 

T2C0F  *  1.5D0  *  Cl 

XLCOF  =  . 125D0  *  A30VK2*SINI0* ( 3 . DO  ♦  5.00  *  COSIO)/(1.DO  +  COSIO) 
AYCOF  *  . 25D0  *  A30VK2  *  SINIO 
X7THM1  *  7. DO  *  THETA2  -  1.00 

90  IFLAG  =  0 


CALL  OPINIT(EOSQ.SINIO, COSIO, BET  AO, A0DP.THETA2, 

1  SING , COSG .BETA02 .XMDOT .OMGDOT .XNODOT .XNODP ) 

C 
C 

C  Update  for  secular  gravity  and 

C  atmospheric  drag. 

C 

c 

100  XMOF  *  XMO  +  XMDOT  *  TSINCE 
C 

OMGADF  *  OMEGAO  +  OMGDOT  *  TSINCE 


C 

XNODDF  «  XNODEO  +  XNODOT  *  TSINCE 
C 


TSQ  «  TSINCE  •  TSINCE 


XNQDE  =  XNODDF  +  XNODCF  *  TSQ 
TEMPA  =  1. DO  -  Cl  *  TSINCE 
TEMPE  *  BSTAR  *  C4  *  TSINCE 
TEMPL  ■  T2C0F  •  TSQ 
XN  ■  XNODP 

CALL  DPSEC( XMDF, OMGADF, XNODE, EM, XINC.XN, TSINCE) 
A  *  (XKE/XN)**TOTHRD  •  T£MPA**2.DO 
E  *  EM  -  TEMPE 

XMAM  *  XMDF  +  XNOOP  *  TEMPL 

CALL  DPPER( E ,XINC .OMGADF .XNQDE , XMAM) 

XL  =>  XMAM  +  OMGADF  ♦  XNQDE 
BETA  »  DSQRT ( 1 .DO  -  E  •  E) 

XN  *  XKE  /  A  ••  1.5DO 

Long  period  penodlcs. 

AXN  1  E  *  DCOS  (OMGADF) 

TEMP  *  1 . DO  /  (  A  •  BETA  *  BETA  ) 

XLL  *  TEMP  *  XLCOF  *  AXN 
AYNL  *  TEMP  *  AYCOF 

XLT  ■  XL  +  XLL 

AYN  *  E  •  DSIN( OMGADF )  +  AYNL 

Solve  Keplers  equation. 

CAPU  *  FMOD2P(  XLT  -  XNODE  ) 

TEMP2  ■  CAPU 
DO  130  I  *  1  , 10 


SINEPV  «  DSIN( TEMP2  ) 


COSEPW  *  DCOS(TEMP2 ) 


TEMP3  =  AXN  *  SINEPW 


TEMP4  *  AYN  *  COSEPW 


TEMP5  «  AXN  *  COSEPW 


TEMP6  *  AYN  *  SINEPW 


EPW  =  ( CAPU  -TEMP4+TEMP3-TEMP2)  /  (1.DO  -  TEMP5  -  TEMP6)* 


I F ( DABS ( EPW  -  TEMP2)  .LE.  ESA )  GO  TO  140 


130  TEMP2  =  EPW 


Short  period  preliminary  quantities. 


140  ECOSE  *  TEMP5  +  TEMP6 


ESINE  =  TEMP3  -  TEMP4 


ELSO  *  AXN  *  AXN  +  AYN  •  AYN 


TEMP  «  1 . DO  ■  ELSO 


PL  =  A  *  TEMP 


R  «  A  •  ( 1 .DO  -  ECOSE ) 


TEMPI  =  1 .DO  /  R 


RDOT  «  XKE  *  OSQRT ( A  i  *  ESINE  *  TEMPI 


RFDOT  «  XKE  *  DSORT(PL)  *  TEMPI 


TEMP2  *  A  *  TEMPI 


BETAL  *  DSQRT (TEMP ) 


TEMP3  «  1.00/(1.00*  BETAL ] 


COSU  *  TEMP2  •  (COSEPW  -  AXN  +  AYN  *  ESINE  *  TEMP3 ) 


SINU  *  TEMP2  *  (SINEPW  -  AYN  -  AXN  *  ESINE  *  TEMP3) 


U  '  ACTAN(  SINU.  COSU  ) 


SIN2U  *  2.00  •  SINU  •  COSU 


C0S2U  «  2.00  •  COSU  *  COSU  -  1.00 


,»  t,»  i.t  ».♦*< 


TEMP  '  1.00/  PL 


TEMPI  *  CK2  *  TEMP 


TEMP2  *  TEMPI  *  TEMP 


Update  for  short  periodlcs  If  required  by  input  card. 


IF  (  IUPSP  .EQ.  'OSCU' )  THEN 


RK  =  R*( 1.D0- 1 .500*TEMP2«B£TAL*X3THM1 )+ . 5D0*TEMP 1 *X 1MTH2 »CQS2U 


UK  =  U  -  .2500  *  TEMP2  *  X7THM1  *  SIN2U 


XNODEK  «  XNODE  +  1.500  *  TEMP2  *  COSID  *  SIN2U 


XINCK  =  XINC  +  1.5D0  *  TEMP2  *  COSIO  *  SINIO  *  CDS2U 


RDOTK  *  ROOT  -  XN  *  TEMPI  *  X1MTH2  *  SIN2U 


RFDOTK  *  RFDOT  +  XN  •  TEMPI  *  ( X1MTH2  *  C0S2U  +  1.5D0  *  X3THM1 ) 


UK  *  U 


XNODEK  *  XNODE 


XINCK  =  XINC 


RDOTK  «  ROOT 


RFDOTK  =  RFDOT 


Orientation  vectors. 


5INUK  *  DSIN(UK) 


COSUK  •  DCOS(UK) 


SINIK  *  DSIN(XINCK) 


COSIK  «  DCOS(XINCK) 


SINNOK  «  DSIN( XNODEK ) 


COSNOK  «  DCOS( XNODEK) 


XMX  «  -  SINNOK  *  COSIK 


XMY  » 


COSNOK  •  COSIK 


C 

C 

C 

C 

C 

r 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


UX  *  XMX  *  S1NUK  +  COSNOK  *  COSUK 

UY  *  XMY  *  SINUK  +  SINNOK  *  COSUK 

UZ  *  SINIK  •  SINUK 
VX  *  XMX  *  COSUK  -  COSNOK  *  SINUK 

V Y  «  XMY  *  COSUK  -  SINNOK  *  SINUK 

V Z  *  SINIK  *  COSUK 

Position  and  velocity. 


X  *  RK  •  UX 

Y  *  RK  *  UY 

Z  *  RK  •  UZ 

XDOT  *  RDOTK  •  UX  +  RFDOTK  *  VX 

YDOT  «  RDOTK  •  UY  +  RFDOTK  •  VY 

ZDOT  «  RDOTK  •  UZ  +  RFDOTK  »  VZ 

RETURN 

END 
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PROGRAM  RDORB1 


FUNCTION 


This  program  reads  a  GTDS  0RB1  file  and  converts  position  and 
velocity  to  Keplerian  and  equinoctial  elements  and  then  writes 
the  elements  into  a  data  set. 


INPUT 


FILE  24  GTDS  epnemeris  run  orbital  elements,  e.g., 

ADC . 0RB1  .DATA 

FILE  5  Gravitational  constant  of  GTDS  EPHEM  ORBl  file 

OUTPUT  . . . . . . . . . . 

FILE  4  Temporary  file,  e.g., 

ADC . EPHEM. DATA 

SUBROUTINES  called  . . . . . 

OPNORB  RDORB  EOUIN  KEPEQN 


HISTORY 


VERSION:  September  1986 

Fortran  subroutine 

ANALYSIS 

Mart i n  E ,  F leger  -  - 

PROGRAMMER 

Mart  in  E .  F ieger  -- 


for  the  IBM  3081  and  3033. 


Captain,  USAF/AFIT/MIT 


Captain,  USAF/AFIT/MIT 


DECLARATIONS 


s  *» 


1  i^iiLk*  !.<  ««*  I.«  J4  AV4  ■»*_«**  y Li ‘ 


IMPLICIT 


LOGICAL 


DOUBLE  PRECISION  (A-H.O-Z) 


Dimensions  »**»*»*»***»»»«**«»*»*«»•***•**»*»»*«••»»••**•»», 


DIMENSION 


P0S( 3 )  ,  VEL( 3 ) 


DOUBLE  PRECISION  EQNELM(6),  KEPELM(6) 


. . DATA  STATEMENTS  »»••*•**«**•*»»***»*»»»»*•»**•«» 


DE2RA  /  .  174532925D- 1  / 


<>.<>>>■>■>>>>>■  BEGIN  PROGRAM  »••»*»»»*»»**»*»*»*•»*»»»***»»*»**» 

Read  the  gravitational  constant  used  on  the 
ORB  1  file. 

READ  (5,900)  GM 

Open  the  GTDS-generated  0RB1  file  24. 

CALL  OPNORB  (  DAYBEG,  SECBEG,  DAYEND,  SECEND,  TIMDIF, 

•  ICENT,  ICOORD.  OAYREP ,  SECREF,  DONE,  24  ) 

Read  another  data  reco’d. 

100  CONTINUE 

CALL  RDORB  (  TIME,  POS,  VEL,  DONE  ) 

IF  (DONE)  STOP 

Calculate  the  Julian  date  of  the  record 
XJUL  =  DAYBEG  +  (  SECBEG  +  TIME)  /  86400. DO 

Calculate  the  YYMMDDHHMMS S .  form  of  the  date. 


(  XJUL  -  IJUL  )  *  86400. DO 


CALL  CALPAK  (YMD.HMS,  DAY JUL , SECJUL ) 


YMD  *  1 .06  ♦  HMS 


Save  the  position  and  velocity  vectors. 


*  POS ( 1  ) 


J-.-.-.v;/.  v.- v„- .-.v-v.v. ■.-.v. A.* 


y  =  pos(2  ) 

2  =  POS ( 3  ) 

OX  =  VEl(1) 

DV  =  VELl 2  ) 

DZ  =  VE  L  ( 3  ) 

Convert  to  equinoctial  elements. 

RETRO  *  1 . DO 

CALL  EQUIN  (  EQNELM,  RETRO,  POS,  VEL,  GM,  .TRUE.  ) 

SMA  *  EQNELM(I) 

XH  *  EQNELM(2 ) 

XK  =  E0NELM( 3 ) 

XP  =  EQNELM( 4 ) 

XO  *  EQNELM( 5 ) 

XML  =  EQNE  LM( 6 ) 

Convert  to  Keplerlan  elements. 

CALL  KEPEQN  (  KEPELM,  EQNELM,  RETRO  ) 

ECC  =  KEPELM( 2 ) 

XINC  =  KEPELM( 3 ) 

XLAN  *  KEPELM( 4 ) 

AP  =  KEPELM( 5 ) 

XMA  =  KEPELM(6) 

Calculate  the  perigee  radius  and  apogee  radius. 

PR  *  SMA  •  (  1.DO  -  ECC  ) 

APR  =  SMA  *  (  1 .00  +  ECC  ) 

Convert  angular  quantities  to  degrees. 


XML 

= 

XML 

/ 

0E2RA 

XINC 

= 

XINC 

/ 

DE2RA 

XLAN 

= 

XLAN 

/ 

DE2RA 

AP 

S 

AP 

/ 

DE2RA 

XMA 

Z 

XMA 

/ 

DE2RA 

Write  tne  Julian  day  and  orbital  elements. 

WRITE  (4,1000)  XJUL,EPOKE,X,V,Z,DX,DY,DZ, SMA , PR , APR , ECC , 
•  XINC , XLAN , AP ,XMA,XH,XK,XP,XQ,XML 


GO  TO  100 

a................  F0rmat  statements 

900  FORMAT ( 3 1 X , F 10 . 2 ) 


1000  FORMAT ( 'DAY ' ,G19. 10.8X ,G21 . 15 
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w  iw  w  iwwww.'wiiwvw  *j»  ■»  ’■»  ■>  »ji  vj  ».v  *>■  *jrv7jrrjrT  v-»v  vt.vi.tvt 


'.wri.vinM 


•  'X  ' ,G19. 10,7X,'Y  ' ,G19. 10,7X,'Z  '.G19.10  / 

*  'DX  ,G19. 10.7X, 'DY  ' ,G19. 10.7X,'0Z  '.G19.10  / 

*  'SMA' ,G19. 10.7X. 'PR  ' ,G19 . 10, 7X , ' APR ' ,G19 . 10  / 

•  'ECC' ,G19.10,7X,'INC' ,G19.10,7X,'LAN' ,G19. 10  / 

■  'AP  ' , G19 . 10.7X , ' MA  ' ,G19. 10,7X,'H  '.G19.10  / 

•  'K  '  ,G19. 10,7X,'P  '  ,G19.10,7X,'0  '.G19.10  / 

*  'ML  ' ,619 . 10  ) 

C 

END 
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.OV'.V 


X' 


,\'.V .s 


...  -  .Y. 


PROGRAM 


PLOTTER 


C 
C 
C 
C 

C  FUNCTION 

C  . . 

c 

c 


C  This  program  reads  a  set  of  nineteen  orbital  elements  from  a 

C  file  of  GTDS  EPHEM  program  elements  and  performs  a  five-point 
C  Lagrangian  interpolation  of  the  elements  to  the  times  of  a  set  of 
C  NORAD  mean  elements  (which  served  as  the  GTDS  DBS  card  deck  I  . 

C  In  a  general  application,  this  program  reads  a  data  set  of  up  to 

C  nineteen  functions  of  a  single  variable,  which  need  not  be  given 
C  at  regular  intervals;  reads  a  second  data  set  containing  nineteen 
C  functions  of  the  same  variable;  and  interpolates  the  values  of  the 

C  functions  of  the  second  data  set  to  the  values  of  the  independent 

C  variable  of  the  first  data  set. 

C 

C  . 

C  USAGE 

C  ..... 

c 

C  This  program  creates  the  formatted  data  set  and  the  control  card 

C  deck  to  perform  the  CSDL  PL0T4B  program. 

C  The  GTDS  EPHEM  interval  is  the  interval  (in  seconds)  between  data 

C  points  on  the  ORBl  file. 

C 

c 

C  SUBROUTINES  CALLED  . . . * . 

C 

C  LAGRAN  INTERP  RCAERR  LINANG  CARTES 

C 

C-** . * .  HISTORY  . . . * . . 

C 

C  VERSION:  OCTOBER  1986 

C  Fortran  program  for  the  IBM  3081  and  3033. 

C 

C  ANALYSIS 

C  Martin  £  Fieger,  Capt,  USAF  --  AFIT  /  MIT 

C 

C  PROGRAMMER 

C  Martin  E  Fieger,  Capt,  USAF  --  AFIT  /  MIT 

C 

C 

. . * .  DECLARATIONS  ••••*•*••••••••••••••••«•»»••••••••*•••• 


C 

c 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C 

C  Character  variables  •***•••**•>•«••**»»•»»« 

C 

CHARACTER  *  7  SVDESG 
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CHARACTER  • 

CHARACTER  * 

CHARACTER  • 

CHARACTER  ■ 

CHARACTER  * 

7 

28 

28 

7 

5 

ytitle 

MTITLE 

DTITLE 

ETITLE 

STITLE 

(20) 

(20) 

(3) 

(3) 

(2) 

Dimensions 

DIMENSION 

DAYCNT 

(5  ) 

DIMENSION 

COE  F 

(11.5). 

. 

(11,5), 

DIMENSION 

POS 

(3). 

VEL 

(3) 

DIMENSION 

ELEM 

(  1000,23,2), 

DIFF 

( 1000,20 : 

DIMENSION 

XMEAN 

(20) , 

SIGMA 

(20) 

DIMENSION 

VECT1 

(6). 

VECT2 

(6). 

DIMENSION 

ANGLIN 

(4,5), 

ANGSAW 

(4,5), 

DIMENSION 

EONELM 

(6) 

c 

C  Data  statements 

C 


DATA 

MTITLE 

/  ' 

1 

X-COMPONENT  OF  POSITION 

1 

Y-COMPONENT  OF  POSITION  ' 

* 

Z-COMPONENT  OF  POSITION  ' 

t 

X-COMPONENT  OF  VELOCITY 

t 

'  Y-COMPONENT  OF  VELOCITY  ' 

• 

Z-COMPONENT  OF  VELOCITY  ' 

• 

'  SEMIMAJOR  AXIS 

» 

'  RADIUS  OF  PERIGEE 

t 

'  RADIUS  OF  APOGEE 

t 

'  ECCENTRICITY 

t 

'  INCLINATION 

* 

'  LONGITUDE  OF  ASCENDING  NODE 

1 

'  ARGUMENT  OF  PERIGEE 

1 

MEAN  ANOMALY 

» 

H 

1 

K 

t 

P 

» 

0 

» 

'  MEAN  LONGITUDE 

/ 

DATA 

YTITLE 

!  -  -  t 

'KM 

'KM 

'KM 

'KM/SEC  '  , 

'KM/SEC  '  , 

KM/SEC  '  , 

'KM 

'KM 

'KM 

POLY  (11) 


ERROR  (6) 
ANGVE  L ( 4 ) 


'DEGREES ' 
'DEGREES' 
'DEGREES 


'DEGREES 


? 


C 


C 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


'DEGREES'  / 


DATA  DTITLE  /  '  POSITION  RADIAL  ERROR  '  , 

'  POSITION  CROSS-TRACK  ERROR  '  , 
'  POSITION  ALONG-TRACK  ERROR  '  / 


DATA  ETITLE  /  'KM  '  , 

'KM 

'KM  '  / 

DATA  TWOPI  /  6.2831  85307  17958  65  DO  / 

DATA  DE2RA  /  0.1745  3292  5  D-1  / 


Begin  program 


INPUT  CONTROL  CARDS 


Read  satellite  designator,  the  GTDS  EPHEM  program 
time  interval,  the  comparison  start  date,  the  first 
day  of  the  plot  and  the  final  day  of  the  plot  (in 
days  since  the  start  date),  the  Fortran  file  numDer 
of  the  irregularly-spaced  data  and  the  legend  name  to 
be  used  on  the  plots,  and  the  Fortran  file  number 
of  the  regularly-spaced  data  and  the  legend  name  to 
be  used  on  the  plots . 

READ  (5,2000)  SVOESG, GINTVL, START, FIRST, FINAL, IRN1.STITLEO), 

IRN2 ,STITLE(2  ) 

Read  the  gravitational  constant  used  on  the 
ORB  1  file. 

READ  (5,2100)  GM 

Convert  the  interval  to  days. 

HAFWID  *  (  GINTVL  *  2 . DO  )  /  86400. DO 

Convert  the  start  date  to  a  Julian  date. 

CALL  JULPAK  (  DAY JUL , SECJUL ,  START ,0. DO  ) 

DAYREF  =  DAY JUL  +  SECJUL  /  86400. DO 
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INTERPOLATOR 

**»*«*****»» 

Set  flags  and  counters  to  initial  conditions. 


Read  in  first  data  record  from  the  GTDS  predictions 
on  file  IRN2. 

READ  ( I RN2  , 3000 , END  =  800 )  DAYCNT ( 1 ) ,  X  ,  YY  , Z  ,DX  ,DY  ,DZ , SMA , PR , APR , 

•  ECC, XINC, XLAN, AP ,  XMA , XH , XK ,XP , XQ , XML 

Save  the  Keplerian  and  equinoctial  elements  for 
interpolation;  convert  angular  quantities  to  radians 


Y  ( 

1,1) 

x 

SMA 

Y( 

2,1) 

= 

ECC 

Y  ( 

3,1  ) 

= 

XINC  »  DE2RA 

Y  ( 

4,1) 

= 

XLAN  •  DE2RA 

Y  ( 

5,1) 

X 

AP  '  DE2RA 

Y  ( 

6,1  ) 

= 

XMA  *  DE2RA 

Y( 

7,1  ) 

x 

XH 

Y  ( 

8,1  ) 

= 

XK 

Y  ( 

9,1) 

= 

XP 

Y  f  10,  1 ) 

£ 

XO 

Y( 11  ,1  ) 

s 

XML  *  DE2RA 

100  CONTINUE 


Read  in  the  next  four  records. 


DO  200  0=2,5 

READ  ( I RN2 , 3000 , END  =  800 )  DAYCNT(J),  X  ,  Y Y  ,  2  ,DX  ,DY ,DZ , SMA ,PR , APR , 

•  ECC, X  INC, XLAN, AP, XMA, XH,XK,XP,XQ, XML 

Save  the  Keplerian  and  equinoctial  elements  for 
interpolation;  convert  angular  quantities  to  radians 


Y<  1,0) 

£ 

SMA 

Y (  2.0) 

£ 

ECC 

Y (  3.0) 

£ 

XINC 

*  0E2RA 

Y (  4,0) 

E 

XLAN 

*  DE2RA 

Y  f  5,0) 

= 

AP 

•  DE2RA 

Y (  6,0) 

£ 

XMA 

*  DE2RA 

Y<  7,0) 

£ 

XH 

Y (  8,0) 

E 

XK 

v<  9.0) 

E 

XP 

Y<  10,0) 

£ 

XO 

Y( 1 1 ,0) 

E 

XML 

*  DE2RA 

yc>:vcyc 


200 

CONTINUE 

c 

c 

Calculate  the  angular  velocity  of  mean  anomaly 

c 

ana  mean  longitude. 

c 

AVGSMA 

=  (  (  Y(1.1)  +  Y( 1  , 5 )  )  /  2. DO  )  *•  3. DO 

ANGVE  L ( 1 ) 

*  DSQRT  (GM  /  AVGSMA ) 

ANGVEU2) 

=  ANGVE  L ( 1  ) 

c 

c 

To  calculate  the  linearized  values  of  the  angular 

c 

elements,  assume  the  longitude  of  ascending  node 

c 

and  argument  of  perigee  are  constant. 

c 

ANGVEL ( 3 ) 

=  0.00 

ANGVEL ( 4 ) 

=  0 .  DO 

c 

c 

Store  the  sawtooth  values  of  the  angular  values. 

c 

DO  250  0*1 .5 

c 

ANGSAWI  1 , 

J)  *  Y(  6,0) 

ANGSAWl 2 , 

0)  =  Y( i 1 ,0 ) 

ANGSAWI 3 i 

0  )  *  Y(  4 ,0  ) 

ANGSAWU, 

0)  *  Y (  5,0) 

c 

250 

CONTINUE 

c 

c 

Calculate  the  linearized  values  of  the  angles. 

c 

CALL  LINANG  ( 

ANGLIN,  ANGSAW,  ANGVEL,  GINTVl,  4,  5) 

c 

c 

Store  the  new  values  of  tne  angular  values  in  the 

c 

array  Y . 

c 

DO  270  0  =  1  ,5 

c 

V(  6,0) 

«  ANGLINf 1,0) 

V  (  1 1 ,  U ) 

*  ANGLIN(2,U) 

v(  4,0) 

*  ANGLINI3.0) 

V(  5.0) 

=  ANGLIN! 4,0) 

c 

270 

CONTINUE 

c 

Set  a  flag  to  snow  interpolator  coefficients  nave 

c 

not  been  calculated 

c 

ICO  EE 

0 

c 

c 

Calculate  tne  center  and  upper  bound  of  tne 

c 

interpolation  interval. 

c 

CENTER 

(  DAYCNT(5)+DAYCNT( 1  )  )  /  2. DO 

c 

UPRBND  = 

CENTER  +  HAFWID 

c 
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c 

c 

c 


If  the  interpolation  interval  is  being  shifted, 
omit  the  next  read  step  (data  was  previously  read). 

IF  (ISHIFT  .EQ.  1 )  GO  TO  400 


C 

C  Read  in  the  NORAD  observations  on  file  IRN1 . 

C 

300  READ  ( I RN 1 ,3000 , END '900 )  DAYADC , X . YY , Z ,DX ,DY ,DZ ,SMA , PR , APR , 

*  ECC , XINC , XIAN , AP , XMA ,H , XK , P ,Q , XML 

C 

C  Check  that  the  date  of  the  data  is  within  the  time 

C  period  for  this  plot. 

C 

IF  (  (DAYADC-DAYREF )  . LT .  FIRST  )  GO  TO  300 
C 

IF  (  (DAYADC-DAYREF)  . GT .  FINAL  )  GC  TO  900 
C 

C  Increment  the  counter. 

C 


C 

C 

C 

C 


c 

c 

c 

c 

c 

c 

c 

c 

c 


NUMADC  =  NUMADC  ♦  1 

Convert  angular  values  to  radians  and  store  all 
the  NORAD  mean  elements  in  array  ELEM. 


ELEM 

( 

NUMADC,  1,1  ) 

= 

DAYADC 

-  dayref 

ELEM 

( 

NUMADC,  2,1  ) 

s 

X 

ELEM 

( 

NUMADC.  3,1  1 

2 

YY 

ELEM 

( 

NUMADC,  4,1  1 

S 

Z 

ELEM 

( 

NUMADC,  5,1  ) 

s 

DX 

ELEM 

( 

NUMADC.  6,1  ) 

= 

DY 

ELEM 

( 

NUMADC,  7,1  ) 

DZ 

ELEM 

( 

NUMADC,  8,1  i 

SMA 

ELEM 

( 

NUMADC,  9,1  1 

= 

PR 

ELEM 

( 

NUMADC, 10,1  > 

=■ 

APR 

ELEM 

( 

NUMADC, 11,1  ) 

* 

ECC 

ELEM 

( 

NUMADC, 12,1  ) 

XINC  * 

DE2RA 

ELEM 

( 

NUMADC, 13,1  ) 

= 

X LAN  • 

DE2RA 

ELEM 

( 

NUMADC, 14,1  ) 

AP  * 

DE2RA 

ELEM 

( 

NUMADC, 15,1  ) 

= 

XMA  * 

DE2RA 

ELEM 

( 

NUMADC, 16,1  ) 

z 

H 

ELEM 

( 

NUMADC, 17,1  ) 

= 

XK 

ELEM 

( 

NUMADC, 18,1  ) 

= 

P 

ELEM 

( 

NUMADC, 19,1  ) 

= 

0 

ELEM 

( 

NUMADC, 20,1  ) 

= 

XML  * 

DE2RA 

Is  the  time  within  bounds  of  the  current 
interpolator  interval"  If  so,  calculate  the 
independent  variable  for  the  interpolator  and 
increment  the  counter.  If  not,  shift  the 
interpolator  interval  by  saving  the  last  element 
set  and  reading  four  more  element  sets.  Set  a  flag 
to  denote  a  shift  is  in  progress. 


400  IF  (  DAYADC  . L E  UPR8ND  )  GO  TO  600 


I  SHI  FT 


1 


DAYCNT  (  1  )  =  DAYCNT(5) 

DO  500  1=1.11 

Y  ( I  ,  1  )  =  Y  ( I  ,  5  ) 

500  CONTINUE 

GO  TO  100 

600  CONTINUE 

A  new  interpolator  interval  Has  been  found. 
Calculate  the  new'  value  of  the  independent 
variable  and  reset  the  shift  flag. 

X VAR  =  ( DAYADC  -  CENTER  )  /  HAFWID 

ISHIFT  =  0 

Calculate  the  interpolator  coefficients  for  this 
period  if  the  coefficient  flag  is  off.  Then  turn 
the  flag  on. 

IF  (  ICOEF  . EO -  0  )  THEN 

CALL  LAGRAN  (  COEF,  Y.  11 .  5  ) 

ICOEF  =  1 

END  IF 

Interpolate  the  equinoctial  anO  Keplerian 
elements  for  the  date  DAYADC. 

CALL  INTERP  (  POLY.  COEF,  XVAR ,  11,  5  ) 

Convert  from  equinoctial  elements  to  Cartesian 
elements . 


EONELM 

(  1  ) 

* 

POLY 

(  1) 

EONELM 

(2) 

S 

POLY 

(  7) 

EONELM 

(3) 

S 

POLY 

(  8) 

EONELM 

(4) 

s 

POLY 

(  9) 

EONELM 

(5) 

= 

POLY 

(  10) 

EONELM 

(6) 

z 

POLY 

(ID 

RETRO 

- 

1.00 

CALL  CARTES  (POS,  VEL,  EONELM,  RETRO,  GM  ) 


Store  the  interpolated  elements  in  array  ELEM 


ELEM 

(NUMADC,  2,2) 

* 

POS( 1 ) 

ELEM 

(NUMADC,  3,21 

= 

POS ( 2 ) 

ELEM 

(NUMADC,  4,2) 

« 

POS( 3  ) 

ELEM 

(NUMADC,  5,2) 

= 

VE  LIU 

ELEM 

(NUMADC,  6,2) 

VEL ( 2  ) 

ELEM 

(NUMADC.  7,2) 

VE  L ( 3  ) 

ELEM 

(NUMADC,  6,2) 

= 

POLY(  1  ) 

ELEM 

(NUMADC,  9,2) 

* 

P0LY(1)  * 

(  1 .DO  -  POL Y ( 2 ) 

ELEM 

(NUMADC, 10.2 ) 

= 

POLY (  1  )  » 

(  1 .DO  +  POLY (2 ) 

ELEM 

(NUMADC, 1 1 ,2 ) 

r 

P0LY(2) 

ELEM 

(NUMADC, 12 .2 ) 

= 

POL  Y ( 3 ) 

ELEM 

(NUMADC, 13,2) 

= 

POLY ( 4  ) 

ELEM 

(NUMADC, 14, 2) 

= 

POL Y( 5  ) 

ELEM 

(NUMADC, 15,2 ) 

= 

POL Y ( 6  ) 

ELEM 

(NUMADC, 16 ,2 ) 

= 

POLY ( 7  ) 

ELEM 

(NUMADC, 17,2 ) 

* 

POLY ( 8 ) 

ELEM 

(NUMADC ,18,2) 

= 

POLY ( 9 ) 

ELEM 

(NUMADC, 19,2  I 

POLY ( 1C) 

ELEM 

(NUMADC, 20, 2) 

POLY ( 1 1 ) 

Calculate  the  differences  Detween  the  GTDS 
predictions  and  the  NORAD  mean  elements. 


DO  660  J=2,20 

DIFF(NUMADC,JI  «  ELEM(NUMADC,J,2  I  -  ELEMI  NUMADC ,J  ,  1  ) 

Filter  tne  angular  elements.... 

IF  (  J.E0.12  .OR  J.EG.13  .OR.  J.EG.14 

OR .  J.EO. 15  .OR .  0. E0.20  )  THEN 

Calculate  the  modulus  by  two  pi  and  convert  to 
degrees 


I  =  ELEMI NUMADC ,0,1)  /  TWOPI 

ELEMINUMAOC.J, 1  )  =  E  LEM ( NUMADC ,0 , 1 )  -  I  •  TWOPI 

E  LEM( NUMADC , U , 1 )  »  ElEM(NUMADC,J,1  )  /  DE2RA 

IF  (  ELEMINUMAOC.J,  1 )  . LT .  0 . DC  ) 

ELEMINUMAOC.J, 1 )  =  E LEM( NUMADC ,J , 1 )  +  360 . DO 

Calculate  tne  modulus  Dy  two  pi  and  convert  to 
degrees . 

I  =  ELEMINUMAOC.J, 2)  /  TWOPI 

ELEMI NUMADC ,J, 2 )  =  E LEM( NUMADC ,J , 2 )  -  I  *  TWOPI 

ELEMINUMAOC.J, 2 )  *  ELEMI NUMADC ,J, 2 )  /  DE2RA 


IF  (  ELEMINUMAOC.J, 2 )  . LT .  0 . DC  ! 


c 

c 

c 

c 


c 

c 

c 

c 

c 


£  L£M( NUMADC ,U  ,2  )  ■=  ELEM(NUMADC  ,  U,2  )  +  360  .  DO 

Calculate  the  modulus  by  two  pi  and  convert  to 
degrees . 

I  =  DIFF(NUMADC.J)  /  TWOPI 

DI FF ( NUMADC  ,  J )  =  DIFF( NUMADC ,J )  -  I  *  TWOPI 

D I F F ( NUMADC , J )  =  DI FF{ NUMADC ,J )  /  DE2RA 


Constrain  the  difference  to  be  between 
minus  pi  and  plus  pi. 


IF 

(  DIFF(NUMADC,J) 

.GT. 

180. DO  ) 

DI FF ( NUMADC ,U ) 

- 

DI FF{ NUMADC , J )  -  360. DO 

IF 

(  DI FF( NUMADC ,J ) 

.LT  . 

-180. DO  ) 

DI FF ( NUMADC , U ) 

DI FF ( NUMADC , U )  +  360. DO 

ENDIF 


660  CONTINUE 


Save  tne  four  position  and  velocity  vectors  in  a 
separate  array . 


DO  700  1=1,6 

VECT2(I)  =  ElEM  (NUMADC,  1+1,  2  ) 

VECTl(I)  =  ELEM  (NUMADC,  1+1,  1  ) 

700  CONTINUE 


Call  the  subroutine  RCAERR  to  calculate  the 
radial,  cross-track,  and  alorq-track  errors  in 
position  . 

CALL  RCAERR  (  VECT1,  VECT2 ,  ERROR  ) 

Store  the  errors  in  the  array  ELEM. 


DC  750  1*1,3 

ELEM  (NUMADC,  1+20,  1)  *  ERROR  (I) 

75C  CONTINUE 


Read  the  next  DAVADC  date. 


GO  TO  300 


If  there  are  insufficient  ORB i  data  records  to 


interpolate  for  every  NORAD  mean  element  recora, 
reduce  tne  number  of  element  sets  plotted. 

800  NUMADC  *  NUMADC  -  1 

900  CONTINUE 

PL0T4B  DATA  CARDS 

Write  tne  array  ELEM  into  tne  data  set 
formatted  by  tne  control  card  for  tne 
plotting  program. 

DO  1200  J  =  2 ,20 

XMEAN(d)  *  0 .  DO 

SIGMA(J)  =  0 . DO 

DO  1000  1*1, NUMADC 

Write  element  comparisons  into  file  41. 

WRITE  (41,4000)  ELEM  (1,1,1),  ELEM  (I.U.I).  ELEM  (I.U.2) 
Calculate  tne  mean  difference. 

XMEAN(J)  =  XMEAN(J)  +  DABS  (  ELEM(I,U,2)  -  ELEM(I,J,1)  ) 

1000  CONTINUE 

XMEAN(J)  =  XMEAN(J)  /  NUMADC 

Calculate  tne  standard  deviation. 

DO  1100  1*1 .NUMADC 
SIGMA(U)  =  SIGMA(U)  + 

*  (  DABS(  ELEM(I,U,2)-ELEM(I,U,1)  )  - 

•  XMEAN(U)  )  •*  2 

1100  CONTINUE 

SIGMA(J)  *  DSQRT  (  SIGMA(U)  /  (NUMADC-1)  ' 

1200  CONTINUE 

Write  the  radial,  cross-track,  and  along-track 
errors  into  file  41. 

DO  1300  U  =  21  ,23 


DO  1300  1*1, NUMADC 


WRITE  (41,4100)  ELEM  (1,1,1),  ELEM  (1,0,1] 


1300  CONTINUE 


Write  the  difference  between  the  GTDS-generated 
elements  and  the  NORAD  mean  elements  into  file  4i . 


DO  1400  J-2 ,20 


DO  1400  I'l.NUMAD C 


WRITE  (41,4100)  ELEM  (1,1,1).  DIFF  (I.J) 


1400  CONTINUE 


PL0T4E  CONTROL  CARDS 


write  the  PL0T4B  control  cards  for  the 
comparison  plots 


DO  1500  U=2 ,20 


WRITE  (40,5000)  NUMADC,  SVDESG.  MTITLE(J),  START, 
STITLE(1).STITLE(2), 

XMEAN(U),  SIGMA(J),  NUMADC 

WRITE  (40,5010)  FIRST,  FINAL,  YTITLE(J) 


1500  CONTINUE 


Write  the  PL0T4B  control  cards  for  the 
error  plots 

DO  1600  <J*1  ,3 

WRITE  (40,5100)  NUMADC,  SVDESG,  DTITLE(U),  START, 
*  STITLE( 1  ).STITtE(2) 

WRITE  (40,5110)  FIRST,  FINAL.  ETITLE(U) 

1600  CONTINUE 

Write  the  PL0T4E  control  cards  for  the 
di  f  f erence  plots 

00  1  TOC  U  =  2 ,20 


WRITE  (40,5200)  NUMADC , SVDE SG ,  MTITLE(J),  START, 
STITLE(2),STITLE( 1) 


WRITE  (40,5210)  FIRST.  FINAL,  YTITLE(J) 
1700  CONTINUE 

Write  final  PL0T4B  control  card. 


WRITE  (40,5400) 


STOP 


c«««. ««.»>>*»*>•»»  Format  statements  »»»«*»**»»•»*•»***»»«»*•“ 

C 

C  Control  cards 

C 

2000  FORMAT  (  49X , A7/49X , F7 . 0/49X , F7 . 0/49X , F7 . 0/49X , F7 .0/ 

*  49X,I2/49X,A5/49X,I2/49X,A5  ) 

C 

2100  FORMAT  (  49X.F10.2  ) 

C 

C  Input  files 


C 

3000  FORMAT  (  3X.G19.10/  6  ( 3X , G1 9 . 10 , 10X ,G 19 . 10 , 10X , G 1 9 . 10/  ) 

*  3X.G19.10  ) 

C 

C  PL0T4B  data  cards 

C 

4000  FORMAT  (  F8 . 3 , 2 ( 3X ,G19 . 10 )  ) 

C 

4100  FORMAT  (  F8.3,3X,G19. 10  ) 

C 

C  PLOT 4B  control  cards 

C 

5000  FORMAT 

*  ( ' *DATA  NUMVAR  3 

*  '  NOREWIND  ',14 

*  '  FMTDATA  1 

*  '(  F8.3,2(3X,G19. 10)  ) 

*  '‘TITLE  MAINTITL  0 

*  A7.A28 

*  '  AXTITL  0 

*  'DAYS  ELAPSED  SINCE  ',F7.0 
■  A5 

*  A5 

*  '  SUBTITL 

*  'MEAN  DIFFERENCE :  '  ,  G 1 1 . 5 , 1 X , ' SI GMA :  '.G11.5.1X, 

*  'AFTER  '  ,14,  '  COMPARISONS' 

5010  FORMAT 

*  (-‘PLOTMOD  YSCALE 

*  '  0.  0. 

•  '  XSCALE 


XLlMIT 

, F7 . 0 . '  ,F7.C 

NOZERO  1 

BADPOINT  1002 

SCATTER  1 

'•PLOT  TITLPLOT 

LEGEND  1 

2  3 

BLKPLOT  200 

A7 

112  3 


5100  FORMAT 

•  (''DATA  NUMVAR  2 

•  '  NOREWIND  ',14 

•  '  FMTDATA  1 

■  '(  F6.3.3X.G19. 10  I 

•  "TITLE  MAINTITL  0 

•  A7.A2S 

•  '  AXTITL  0 

•  'DAYS  ELAPSED  SINCE  ,F7.0 
»  ERROR 

•  '  SUBTITL 


*  A5 ,  (ESTIMATED);  ,A5,'  (TRUTH) 


5110  FORMAT 

•  (  "PLOTMOD  YSCALE 

•  '  0.  0- 

*  XSCALE 

•  0-  0. 

•  '  XLIMIT 

•  . F7 . 0 ,  '  ,  F7  .  C 

•  NOZERO  1 

•  '  BADPOINT  1002 

•  SCATTER  1 

•  "PLOT  TITLPLOT 

•  '  LEGEND  1 

•  2 

•  '  BLKPLOT  200 

•  A  7 

•'112 


5200  FORMAT 

•  ("DATA  NUMVAR  2 

•  '  NOREWIND  ',14 

•  '  FMTDATA  1 

•  '(  F8 .3,3X,G19. 10  ) 

•  "TITLE  MAINTITL  0 

•  A7.A28 

•  '  AXTITL  0 

•  'DAYS  ELAPSED  SINCE  '  ,F7 .0 

•  'DELTA 


SUBTITL 


' COMPARISON  DIFFERENCE:  DELTA 


' .65, ' - ' .65 


5210  FORMAT 

■  ( ' 'PLOTMOD  YSC6LE 

*  '  0.  0. 

*  '  XSCALE 

*  '  0.  0. 

*  '  XLIMIT 

*  '  , F7 . 0 , '  '  ,  F7 . 0 

*  '  NOZERO  1 

*  '  BADPOINT  1002 

*  '  SCATTER  1 

*  "PLOT  TITLPLOT 

*  '  LEGEND  1 

*  '  2 

»  BLKPLOT  200 

»  A  7 

*  '  1  1  2 


5400  FORMAT  ("END  LAST  ) 

END 


PROGRAM  RUNADCOB 


P 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FUNCTION 


**»»***» 


This  program  reads  an  ADCOM  observation  card  file  and  writes  a 
GTDS  card  file. 


INPUT 


File  5  ADCOM  observation  cards. 


File  8  Station  acronyms. 

LWE 1122. ORB  IT. ST  AT  FILE. DATA 


OUTPUT  . . 

File  6  Printer  messages. 

File  7  GTOS  observation  cards. 


/OBSDAT / 


c 

Observation 

time 

computed  in  subroutine  DATE. 

c 

c 

DAVOBS 

I 

Uul lan  date  at  noon  on  day  of  observation 

c 

SECOBS 

I 

Time  of  observation  in  seconds  from  noon. 

c 

(Range:  -43200  to  almost  43200) 

c 

c 

YMDOBS 

I 

Calendar  date  packed  in  the  form  YYMMDD 

c 

HMSOBS 

I 

Time  of  day  packed  in  the  form  HHMMSS 

c 

c 

YEAR 

I 

Year  -  1900. 

c 

MONTH 

I 

Month . 

c 

DAY 

I 

Day  . 

c 

c 

HOUR 

I 

Hour  . 

c 

MINUTE 

I 

Minute . 

c 

SECOND 

I 

Second . 

SUBROUTINES  CALLED 
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■S 


a mmmsa 


aWSi*vCii 


i  » 


I'i^l  (.1  tat 


RANGES  A2IMUT  ELEVAT  RANGER  ASTRON 


•»»**»**»**»*  HISTORY 


VERSION:  October  1986 

Fortran  main  program  for  the  IBM  3090. 


ANALYSIS 


Joe  F.  Lombardo 


PROGRAMMER 


Joe  F.  Lombardo 


--  Charles  Stark  Draper  Laboratory 


--  Charles  Stark  Draper  Laboratory 


MODIFIED 
Oct  1986 


Leo  W.  Early,  Jr. 

Charles  Stark  Draper  Laboratory 


1.  Converted  from  FORTRAN  66  to  FORTRAN  77. 

2.  Replaced  algorithm  used  to  find  station  acronym.  Old 
algorithm  searched  entire  station  array.  New  algorithm 
uses  a  single  array  reference. 

3.  Improved  code  structure. 

4.  NEW  --  Print  lists  of  known  stations  used  and  unknown 
stations  used 


MODIFIED 
Dec  1986 


Leo  W.  Early,  Jr. 

Charles  Stark  Draper  Laboratory 


1.  Output  GTDS  observation  cards  in  standard  GTDS  observation 
card  order. 

2.  Transform  coordinates  of  optical  observations  (right  ascen¬ 
sion  and  declination)  from  NORAD  to  mean-of -  1950 . 0  coordinate 
system . 

3.  Suppress  output  of  observation  cards  for  unknown  stations. 


MODIFIED 
Dec  1986 


Leo  W.  Early,  Jr. 

Charles  Stark  Draper  Laboratory 


1.  Add  coordinate  system  for  optical  observations.  See  subrou¬ 
tine  ASTRON. 


DECLARATIONS 


implicit 


double  PRECISION  (A-H.O-ZI 


55W 


/OBSDAT/  »..*«>*.********* 
INCLUDE  ( OBSDAT M 


Constants  •*«**»*«»*»**»*•**»*** 


****»**•***»**»»*«*■********* 


Maximum  station  ID  number. 
PARAMETER  (NSTAT  =  999 


) 


Character 


Variables  »*«»*»•»»»**« 


»»»****»*#**•**** 


CHARACTER  * 

1 

LEL 

■  LRA 

CHARACTER  * 

1 

LELR 

,  LAXR 

CHARACTER  * 

1 

LAC 

,  EQNYR 

character  • 

4 

STANAM 

, STNAME  ( 0 : NST  AT ) 

CHARACTER  ■ 

4 

nulsta 

Numeric  Data  Types  «*•*»»'»«*"**' 
LOGICAL  OBSERV  ( 0 : NST AT  ) 


DATA  STATEMENTS  * 


*#***»»******»********%*•** 


Null  station  acronym. 


DATA 

NULSTA 

/  »****'  j 

Number  of  observations  ignored 

DATA 

NACC 

/  0  / 

DATA 

NCOS 

/  0  / 

DATA 

NCOSR 

/  0  / 

BEGIN  PROGRAM 


READ  STATION  LIST 
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C  Fill  station  list  with  null  acronyms. 

C 

DO  10  ISTNUM  =  O.NSTAT 

STNAME  (ISTNUM)  =  NULSTA 
10  CONTINUE 
C 

C  Read  station  acronyms  from  station  file. 

C 

20  CONTINUE 

READ  (8, 4000, END* 30)  I STNUM . STANAM 

STNAME  (ISTNUM)  =  STANAM 

GO  TO  20 
C 

C  Set  station  oDservation  indicators  --  there  are  no 

C  observations  for  any  station. 

C 

30  DO  40  ISTNUM  =  O.NSTAT 

OBSERV  (ISTNUM)  *  .FALSE  . 

40  CONTINUE 

C 


C 

C 

c 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


* 

k 

N 

I 

ft 


k 

K 

l 


I 

» 

ft 

I 

i 

► 

* 

* 

t 

t 

( 


t, 


c 

c 


Compute  Julian  date  and  GTDS  calendar  date. 


c 

CALL  DATE  ( I YR , IDAV , I HR , IMIN , SEC ) 

C 

C  Compute  GTDS  observation  card  aate. 

C 

OBSTIM  =  YMDOBS  *  1.D6  +  HMSOBS 

C 

C  Find  station  name. 

C 

STANAM  =  STNAME  (ISTNUM) 

C 

C  Branch  to  code  for  this  observation  type. 

C 

GOTO  1200.205.210.215.220,  225.230.235,240,245), 

*  ITYPE  +  1 


C 

C 

C 

C  . . .  WRITE  GTDS  CARDS 

C 

C  .  ACCOM  type  0 . 

C 

C  Range  rata.  (GTDS  type  9) 

C 


200  IF  (STANAM  .NE.  NULSTA)  THEN 

CALL  RANGER  (ILRA.RATR,  RATE,  STANAM ,OBST IM ) 

END  IF 

OBSERV  (ISTNUM)  *  .TRUE. 

GO  TO  100 


C 

C  .  ACCOM  type  1  . 

C 

C  Azimuth.  (GTDS  type  4) 

C  Elevation.  (GTDS  type  5) 

C 


205  IF  (STANAM  .NE.  NULSTA)  THEN 

CALL  AZIMUT  (IAZORA,  AZMTH ,  ST ANAM , OBST IM ) 

CALL  ELEVAT  (ILEL.ELEVR,  ELEV,  ST ANAM , OBST IM ) 

END  IF 

OBSERV  (ISTNUM)  =  .TRUE  . 

GO  TO  100 
C 

C  . ACCOM  type  2 . 

C 


C 

Range . 

(GTDS  type 

1  ) 

C 

Azimuth. 

(GTDS  type 

4  ) 

c 

r 

Elevat ion . 

(GTDS  type 

5) 

L. 

210 

IF  (STANAM  .NE. 

NULSTA)  THEN 

CALL 

RANGES 

( RNGE , I RX  , 

RANGE  , 

STANAM, OBSTIM) 

CALL 

AZIMUT 

(IAZORA, 

AZMTH, 

STANAM, OBSTIM) 

CALL 

ELEVAT 

(ILEL.ELEVR, 

ELEV, 

STANAM, OBSTIM) 

END  IF 

OBSERV  (ISTNUM)  =  .TRUE. 
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GO  TO  100 


R: 


t 


ADCOM  type  3 


215 


220 


225 


Range . 

(GTDS  type 

i  > 

Azimuth. 

(GTDS  type 

4) 

Elevaticn. 

(GTDS  type 

5) 

Range  rate. 

(GTDS  type 

9) 

IF  (STANAM  .NE.  NULSTA)  THEN 

CALL 

RANGES  (RNGE.IRX, 

RANGE . 

STANAM, 0BST1M) 

CALL 

AZIMUT  (IAZORA, 

AZMTH , 

STANAM,  OBSTIM) 

CALL 

ELEVAT  (ILEL.ELEVR, 

ELEV, 

STANAM, OBSTIM) 

CALL 

RANGER  (ILRA.RATR, 

RATE, 

STANAM, OBSTIM) 

END  IF 

OBSERV 

(ISTNUM)  =  .TRUE. 

TO  100 

-  ADCOM  type  A  - 

-- 

Range . 

(GTDS  type 

1  ) 

Azimuth. 

(GTDS  type 

4  ) 

E  1  evat  ion . 

(GTDS  type 

5) 

Range  rate. 

(GTDS  type 

9) 

Range  acceleration. 

Azimuth  rate. 

Elevation  rate. 

IF  (STANAM  .NE.  NULSTA)  THEN 

CALL 

RANGES  (RNGE.IRX, 

RANGE  , 

STANAM, OBSTIM) 

CALL 

AZIMUT  (IAZORA, 

AZMTH, 

STANAM, OBSTIM) 

CALL 

ELEVAT  (ILEL.ELEVR, 

ELEV, 

STANAM, DESTIM) 

CALL 

RANGER  (ILF'-.RATR, 

RATE  , 

STANAM, OBSTIM) 

END  Ic 

OBSERV 

(ISTNUM)  =  .TRUE 

NACC 

=  NACC  +  1 

TO  100 

.  ADCOM  type  5  --- 

... 

Right  ascension. 

(GTDS  type 

6) 

Decl i nat i on . 

(GTDS  type 

7) 

IF  (STANAM  .NE.  NULSTA)  THEN 

CALL 

ASTRON  (IAZORA, ILEL, 

ELEVR,  EONYR,  STANAM) 

END  IF 

OBSERV 

(ISTNUM)  =  .TRUE. 

TO  100 

-  ADCOM  type  6  - 

... 

Range . 

(GTDS  type 

1) 

230 


IF  (STANAM  .NE.  NULSTt)  THEN 
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I 


./W’ 

VJV. 


CALL  RANGES  (RNGE.IRX,  RANGE,  STANAM  .OBSTIM ) 


OBSERV  (ISTNUM) 


GO  TO  100 


ADCOM  type  7 


Direct i on  cosines . 


235  NCOS 
GO  TC  100 


-  ADCOM  type  8 

Range . 

Direction  cosines. 


{ GTDS  type  1) 


240  IF  (STANAM  . NE .  NULSTA )  THEN 

CALL  RANGES  (RNGE.IRX,  RANGE,  STANAM  .OBSTIM ) 


OBSERV  (ISTNUM) 


NCOSR  +  1 


GO  TO  100 


ADCOM  type  9 


Range . 

Range  rate. 

D 1  recti  on  cos i nes . 


(GTDS  type  1) 
(GTDS  type  9) 


245  IF  (STANAM  .NE.  NULSTA)  THEN 

CALL  RANGES  (RNGE.IRX,  RANGE,  STANAM  .DBSTIM ) 

CALL  RANGER  (ILRA.RATR,  RATE,  ST ANAM , OBST IM ) 

END  IF 


OBSERV  (ISTNUM) 


GO  TO  100 


NCOSR  *  1 


END  FILE 


Write  END  card. 


300  WRITE  (7,3000)  'END 


PRINT  STATISTICS 


c 


*  *  *  *  * 


LIST  Of  KNOWN  STATIONS 


C 
C 

C  Print  heading  --  List  of  known  stations  used. 

C 

WRITE  (6,2000) 

LINE  =  5 

C 

C  Are  tnere  observations  for  this  station? 

C  Is  this  station  known? 

C 

DO  410  ISTNUM  =  O.NSTAT 

IF  ( OBSERV  (ISTNUM)  .AND. 

*  STNAME  (ISTNUM)  . NE .  NULSTA)  THEN 

C 

C  Print  blank  line  after  every  five  station  names 

C 

LINE  =  LINE  +  1 
IF  (LINE  .GT.  5)  THEN 

LINE  =  1 

WRITE  (6,2110) 

END  IF 
C 

C  Print  station  number  and  acronym. 

C 

WRITE  (6,2100)  ISTNUM,  STNAME  (ISTNUM) 

END  IF 

410  CONTINUE 


C 

c 

c 

C  * .  LIST  of  UNKNOWN  STATIONS  . 

c 

C  Print  heading  --  List  of  unknown  stations  usee 

C 

WRITE  (6,2010) 

LINE  *  E 


C  A^e  mere  on 

C  Is  IMS  si»* 

C 

DC  420  ISTNUM  =  0 ,NS'A' 

IF  (OESEP.  )  I S T N JM  i 
•  STNAME  I  I 5TNJM  I 

0  Pr-r.t  a-lfc 

LINE  5  l.IN£  .  • 

I'  ‘  *  I  NE  S'  ?  ■ 

l  I  Nf  « 

wc ;  T  i  i 
ENL  If 

0  . .  s  r  a  • 


e-vat  -o-s  f-s  stat  -o-' 

Or  vnk  no»r ' 

an: 

E  0  Nu^  S'  A  i  T~ik, 

me  a f  *  er  eve  -  .  •  ,  vp  s  a '  '  mmes 

’»E  N 

r  -TOP'  ar-C  8  "  '  I  '  .  1 


WRITE  (6,2100)  ISTNUM,  NULSTA 


CONTINUE 


IGNORED  OBSERVATIONS 


Print  number  of  non-GTDS  observations  Ignored. 


WRITE  (6,2200)  NACC .NCOSR .NCOS 


FORMAT  STATEMENTS 


FORMAT ( 


AOCOM  observation  cards. 

T7 ,  13.  12, 13, 12, 12. F5. 3. 

A1.F5.4.1X,  17, IX,  F7.5.I1.1X,  A1.F6.5.1X, 

A 1  ,  F4 . 4  ,  IX  ,  A  1  ,  F4 . 4 , 1 X  ,  A1.F4.4.2X,  II, A1) 

Printer  messages. 


FORMAT( 


Known  Stations  Used  *•***', 


FORMAT ( 


FORMAT ( 
FORMAT ( 
FORMAT ( 


. .  II 

'  This  list  includes  all  known  stations  in  the  ', 

ADCQM  observation  file  which  generated  GTDS-usable'  / 

'  observations.'  ) 

III  '  . .  Unknown  Stations  Used  ****•', 

'*•*»»'  // 

'  This  list  Includes  all  unknown  stations  in  the 
'ADCOM  observation  file  which  generated  GTDS-usable'  / 

'  observations.'  ) 
lx, 15,  10X.A4  ) 

IX  ) 

III  '  **•»**•»**  Non-GTDS  Observations  Found  ****•', 
'***»*'  // 

IX, 15,'  Elevation  rate,  azimuth  rate,  and  range  ', 
'acceleration  triples  were  found.'  / 

IX, 15,'  Direction  cosine  pairs  were  found  on  range/', 
'range-rate/dlrectlon-cosine  cards.'  / 

IX, 15,'  Direction  cosine  only  cards  were  found.'  II 
'  All  non-GTDS  Observations  found  were  ignored.'  ) 


FORMA T( 


GTDS  observation  cards. 


Station  acronyms . 


SUBROUTINE  ASTRON  ( IAZORA , ILEL .ELEVR , 


EONYR,  STANAM) 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FUNCTION 


This  subroutine  writes  an  ACCOM  optical  observation  pair 
(right  ascension  and  declination)  onto  a  pair  of  standard 
GTDS  observation  cards. 


CALLING  SEQUENCE  ...*»•*»***•«»•**•*•*****•**•***'»•***•*•* 
CALL  ASTRON  ( I A20RA  ,  I LE L  ,  ELEVR  ,  EQNYR,  STANAM) 


PARAMETERS 


IAZORA  I 
ILEL  1 
ELEVR  I 


Right  ascension.  (HHMMSS.S) 
Declination  --  first  digit  with  overpunch. 
Declination  --  later  digits.  (D.DDDD) 


EQNYR  I  Coordinate  system  of  observation. 

'0'  NORAD  true  of  date  at  time  of  obser- 
vat ion 


'  1 '  mean  of  1900 . 0 

’  2  •  mean  of  1920 . 0 

3  mean  of  1950.0 

'  4  mean  of  1975 . 0 

'5'  mean  of  2000.0 

'6'  mean  of  1850.0 

'T  mean  of  1805.0 

' 8 '  mean  of  1875 . 0 

'9  mean  of  1960.0 

'X'  NORAD  true  of  date  at  beginning  of 
Bessel lan  year  of  observation 

>Y'  NORAD  true  of  date  at  midnight,  Uan  0, 
of  year  of  observation 


STANAM  I  Station  name. 


/OBSDAT/ 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Observation 

t  ime 

computed  in  subroutine  DATE. 

DAYOBS 

I 

Julian  date  at  noon  on  day  of  observation. 

SECOBS 

I 

Time  of  observation  In  seconds  from  noon. 

(Range:  -43200  to  almost  43200) 

YMDOBS 

I 

Calendar  date  packed  In  the  form  YYMMDD. 

HMSOBS 

I 

Time  of  day  packed  in  the  form  HHMMSS.SSSSS 

YEAR 

I 

Year  -  1900. 

MONTH 

I 

Month. 

DAY 

I 

Day . 

HOUR 

I 

Hour  . 

MINUTE 

I 

Minute . 

SECOND 

I 

Second . 

subroutines  called  »* 

PRENUT  VEMA33 


HISTORY 


VERSION:  December  1986 

Fortran  subroutine  for  the  IBM  3090. 

ANALYSIS 

Leo  w.  Early,  Jr.  --  Charles  Stark  Draper  Laboratory 

PROGRAMMER 

Leo  W.  Early,  Jr .  --  Charles  Stark  Draper  Laboratory 


DECLARATIONS 


IMPLICIT  DOUDLE  PRECISION  (A-H.O-Z) 

/OBSDAT/  . . . . . 

INC1 JOE  (OBSDAT*) 


Constants 


196 


•  .It' 


_<*  H  I>il4  1 


t  l.|  M.’ J*» 


■  w«.t  1<.  a*.,  it.  4U*iL*iLkaL.,tl_ 


PARAMETER 


Number  of  mean-of -  reference  coordinate  systems  used 
for  optical  observations. 

(NC00R0  *  9 


Character  variables 


CHARACTER  •  1  EONYR 

CHARACTER  ■  4  STANAM 


Numeric  Data  Types  **• 
LOGICAL  INIT 

Dimensions  *•**»*»»*»« 


DIMENSION 

DIMENSION 

Saved  Variables 


UNIT  (3) 
ROTATE  (3.31 


, YRCORD  ( NCOORD I 


DATA  STATEMENTS 


. .  MATHEMATICAL  CONSTANTS  . . 

Ratio  of  circumference  of  circle  to  diameter. 

DATA  HALFPI  /  1.5707  96326  79489  66  D  0  / 

DATA  TWOPI  /  6.2831  85307  17958  65  D  0  / 

Degrees  to  radians. 

DATA  DEGREE  /  1.7453  29251  99432  96  D  -2  / 


TIME  PARAMETERS 


Reference  epoch  and  Julian  date  of  mean-of - 1950 . 0 
coordinate  system. 


DATA 

YRREF 

/ 

1950. 

D 

0 

/ 

DATA 

OAYREF 

/ 

24332  82  423 

D 

0 

/ 

DATA 

SECREF 

/ 

0. 

D 

0 

/ 

C  Reference  epochs  of  observation  coordinate  systems 
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< fc:  «>. 


■  >L  j». 


.  iU'iL  iL  iV'il,  «|  «a»  .I  ,  .1.  iL  il. 


DATA 

YRCORD 

(1) 

/ 

1900. 

D 

0 

/ 

DATA 

YRCORD 

(2) 

/ 

1920. 

D 

0 

/ 

DATA 

YRCORD 

(3) 

/ 

1950. 

D 

0 

/ 

DATA 

YRCORD 

(4) 

j 

1975, 

D 

0 

/ 

DATA 

YRCORD 

(5) 

/ 

2000. 

D 

0 

/ 

DATA 

YRCORD 

(6) 

/ 

1850. 

D 

0 

/ 

DATA 

YRCORD 

(7) 

/ 

1855  . 

D 

0 

/ 

DATA 

YRCORD 

(8) 

/ 

1875  . 

D 

0 

/ 

DATA 

YRCORD 

(9) 

/ 

1960. 

D 

0 

/ 

Length 

of 

the  tropical 

yea  r 

DATA 

TROPIC 

/ 

31 

556  925 

.97  47 

D 

DATA  BESACC 


Secular  acceleration  constant  useo  to  convert  from 
Bessel ian  years  to  tropical  years  (inverse  years) 

/  2.345  D  -11  / 


DATA  INIT 


. .  CONTROL  parameters  .......... 

Initialize  tnis  subroutine  on  the  next  call 
/  TRUE.  / 


BEGIN  PROGRAM 


INITIALIZE 


Initialize  subroutine! 

IF  (INIT)  THEN 

Compute  unit  used  to  measure  right  ascension 

TENSEC  •  DEGREE  /  2400.00 

Compute  length  of  tropica'  year  in  ephemeris  oavs 

TROPIC  •  TROPIC  /  86400.00 

Compute  reference  epochs  of  cooromate  systems  in 
years  since  1900 


.%  *.  A'  s  . 


YRREF  *  YRREF  -  1900. DO 

DO  10  1  *  1.NC00RD 

YRCORD  (I)  *  YRCORD  (I)  -  1900. DO 

10  CONTINUE 
C 

C  This  subroutine  has  been  initialized. 

C 

INI T  •  .FALSE 
END  IF 
C 
C 
C 

C  .............. 

C  CONVERT  ANGLES 

C  . 

c 

c 

c 


c  ..........  RIGHT  ascension 

c 

C  Unpack  r ight . 

C 

ISEC  *  IAZORA 

IMIN  *  ISEC  /  100C 

I HOUR  *  IMIN  /  100 

C 

c  Unpack  left 

RAHOUR  «  I  HOUR 

RAMIN  .  IMIN  -  100  ■  IHOUR 

RASEC  *  ISEC  -  100C  •  IMIN 

C 

C  Convert  to  tenth- seconds  of  time 


C 

ALPHA  .  RAHOUR  •  36000-00  *  RAMIN  1  60C . DO  *  RASEC 

C 

C  Convert  to  radians. 

C 

ALPHA  »  ALPHA  •  TENSEC 
C 
C 

c 

c  . .  DECLINATION  .......... 

C 

C  Convert  initial  blank  to  +0 

C 

If  ( ILEL  .  EO  64  1  ILEL  *  240 

C 

C  Convert  initial  minus  sign  to  -  0 

C 

IF  (I LEl  IQ  96)  ILEL  *  208 

C 

C  Decllnatior  is  positive  --  first  digit  numeric 
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IF  (ILE1-  .GE.  240)  THEN 

DELTA  *  (ILEL  -  240)  *  10  +  ELEVR 

Declination  is  negative  --  first  digit  overpunched 


ELSE 

DELTA  *  -  (ILEL  -  208)  *  10  -  ELEVR 

END  IF 


Convert  to  radians. 

DELTA  *  DELTA  *  DEGREE 


COMPUTE  ROTATION  MATRIX 


.  NORAD  TRUE-OF-DATE  .......... 

Does  observation  use  NORAD  true-of-date  coordinate 
system? 

IF  (EONVP*  .  EQ  'O')  THEN 

Compute  time  of  observation  in  ephemeris  days  since 
195C.0 

WARNING  Tms  computation  ignores  leap  seconds. 

TIME  '  (DAYOBS  -  DAYREH  ♦ 

(SECOBS  -  SECREF)  /  86400. DO 

Compute  transformation  matri.  from  mean-of - 1950 . 0 
coordinates  to  NORAD  true-of-date  coordinates. 

CALL  PRENUT  (ROTATE,  TIME FULL ' ) 


.  NORAD  TRUE-OF-BESSELIAN  .......... 

Does  ooservation  use  the  NORAD  true-of-date  coord¬ 
inate  system  at  the  beginning  of  the  Bessel lan 
year  7 

ELSE  IF  ( EON! R  EQ  x '  )  THEN 


Compute  time  at  beginning  of  Besre  an  year  in 
ephemenis  days  since  1950.0 


AVRYR  =  1 . DO  -  BESACC  *  (YRREF  +  YEAR) 

TIME  *  (YEAR  -  YRREF)  •  AVRYR  ■  TROPIC 

Compute  transformation  matrix  from  mean-of - 1950 . 0 
coordinates  to  NORAD  true-of -Bessel lan  coordinates. 

CALL  PRENUT  (ROTATE,  TIME, 'FULL') 


..........  NORAO  TRUE -OF -UAN-0 . 0  .......... 

Does  observation  use  the  NORAD  true-of-date  coord¬ 
inate  system  at  midnight  on  Uan  0  of  the  year  of 
observat 1  on? 

ELSE  IF  ( EONYR  .EQ.  'Y'  )  THEN 

Compute  Uulian  date  of  uan  0.0  of  year  of 
observat ion . 

CALL  UULIAN  ( DAI JUL , SECUUL ,  YEAR , 1 . DO ,C . DO , 

0.00,0.00,0.00) 

Compute  time  of  Uan  0.0  of  year  of  observation  in 
ephemeris  days  since  1950.0. 

WARNING:  This  computation  ignores  leap  seconds 

TIME  *  (DAYUUL  -  DAYREF)  + 

(SECUUL  -  SECREF)  /  86400 . DO 

Compute  transf ormat ion  matrix  from  mean-of - 1950 . 0 
coordinates  to  NORAD  t rue  -  of -uan-0 . 0  coordinates. 

CALL  PRENUT  (ROTATE,  TIME, 'FULL) 


* . .  MEAN-OF-REFERENCE  . . 

Does  observation  use  a  mean-of -ref erence  coordinate 
system? 

ELSE  IF  (EONYR  . GE .  V  .AND.  EONYR  . LE .  '9')  THEN 

Find  reference  epoch  of  observation  coordinate 
system . 

READ  ( EONYR,  (II)'  )  IYEAR 

If  reference  epoch  of  observation  coordinate  system 
is  1950.0  then  leave  observation  unchanged. 


IF  ( YRCORD  (IYEAR)  . £Q  YRREF)  GO  TO  400 

Compute  time  of  reference  epoch  in  ephemeris  days 
since  1950.0. 

AVRYR  *  1 .DO  -  BESACC  *  (YRREF  +  YRCORD  (IYEAR)) 

TIME  »  (YRCORD  (IYEAR)  -  YRREF)  »  AVRYR  •  TROPIC 

Compute  transformation  matrix  from  mean-of - 1950 • 0 
coordinates  to  mean-of -ref erence  coordinates. 

CALL  PRENUT  (ROTATE,  TIME , ' PRECES ' ) 


mo  COORDINATE  SYSTEM 


Write  error  message. 


ELSE 

WRITE  (6,1000)  ALPHA  .YMDOBS ,  DELTA  .HMSOBS , 


Leave  observation  unchanged. 


GO  TO  400 


TRANSFORM  OBSERVATION 


FIND  UNIT  VECTOR 


Compute  cosine  of  declination. 


COS  (DELTA) 


Compute  unit  vector  from  right  ascension  and 
decl  mat  ion . 


UNIT  (1)  *  COS  (ALPHA)  •  COSDEL 
UNIT  (2)  «  SIN  (ALPHA)  *  COSDEL 
UNIT  (3)  *  SIN  (DELTA) 


TRANSFORM  COORDINATES 


Transform  unit  vector  from  observation  coordinate 


.  •  .  i  T  O 

•  .  *  .*.  *  k  "  »  *  * '  » 


V*i 


«o: 


system  to  mean-of - 1950.0  coordinate  system 
CALL  VEMA33  (UNIT.  UNIT, ROTATE) 


. .  FIND  ANGLES  **••***•*■ 

Compute  magnituoe  of  eouatonal  component  of  unit 
vector . 

UNIEOU  *  UNIT  (1)  *  UNIT  (1)  *  UNIT  (2)  *  UNIT  (2) 

UNI EQU  =  SORT  (UNIEOU) 

Compute  right  ascension  and  declination. 

IF  (UNIEOU  . EO .  O.DO)  THEN 
ALPHA  =  C.OO 

DELTA  =  SIGN  (HALFPI,  UNIT  (3)) 

ELSE 

ALPHA  =  ATANC  (UNIT  (2),  UNIT  (1)) 

IF  (ALPHA  . LT .  C.DOl  ALPHA  =  ALPHA  *  TWOPI 

delta  =■  atan:  ( unit  <3),  unieou; 

END  lc 


WRITE  OBSERVATION  cards 


Compute  GTDS  observation  care  date 
OESTIV  *  VMDOES  •  1.06  *  HMSOES 

write  right  ascension  caro 

WRITE  (T,  20001  STANAM.6,  OESTIM,  ALPHA, ALPHA 

write  declination  cara 

WRITE  (1,20001  STANAM.T.  OBSTIK,  delta, DELTA 

Return 

RETURN 


FORMAT  STATEMENTS 
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c 

c 

c 


Error  message 


1000  FORMAT (  /  '  »*»•***•*•  ERROR  ..........  ODServa 

*  'tton  cara  uses  unknown  coordinate  system  /, 

*  '  Riant  ascension  =  '  ,  1PG17  .  10,20>. ,  Calendar  date 

»  OPFIO.O  / 

*  '  Declination  =  ' , 1PG17 . 10.20X ,  Time  of  da> 

*  OPF13.3  / 

*  '  Coordinate  system  indicator  =  -  ,A  ) 

C 

C  GTDS  observation  cards. 

C 

2000  FORMAT (  A4,5X,  12,6 X,  3G21.14  ) 

END 
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SUBROUTINE  PRENUT  (ROTATE,  TIME, TRANS) 


This  subroutine  computes  precession  and  nutation  matrices  for  the 
Earth  at  a  given  time  and  returns  in  parameter  ROTATE  the 

coordinate  transformation  matrix  requested  in  parameter  TRANS: 


Value  of 
TRANS 


Coordinate  Transformation  in  ROTATE 


r 

'FULL' 

MEAN 

OF 

1950.0 

to 

TRUE 

OF 

DATE 

c 

'PRECES 

MEAN 

OF 

1950.0 

to 

MEAN 

OF 

DATE 

c 

'NUTATE 

MEAN 

OF 

DATE 

to 

TRUE 

OF 

DATE 

CALLING  SEQUENCE 


CALL  PRENUT  (ROTATE,  TIME, TRANS) 


PARAMETERS 


ROTATE  0 


TRANS  I 


Coordinate  transformation  matrix. 

Time  since  1950.0  in  ephemeris  days. 
Coordinate  transformation  requested. 


C 

'FULL' 

MEAN 

OF 

1950.0 

to 

TRUE 

OF 

DATE 

c 

'PRECES  ' 

MEAN 

OF 

1950.0 

to 

MEAN 

OF 

DATE 

c 

'NUTATE 

MEAN 

OF 

DATE 

to 

TRUE 

OF 

DATE 

SUBROUTINES  CALLED 


MSET33  MAMA33 


HISTORY 


VERSION:  February  1985 

Fortran  subroutine  for  the  IBM  3090. 

ANALYSIS 

( unknown ) 


WPS? 


M>vvwV, 


PROGRAMMER 


( unknown ) 


MODIFIED 
Feb  1985 


Leo  W .  Early  ,  or . 

Charles  Stark  Draper  Laboratory 


1.  Converted  code  to  Fortran  77. 

2.  Improved  code  structure. 


DECLARATIONS 


IMPLICIT 


DOUBLE  PRECISION  (A-H.O-Z) 


Character  Variables 


CHARACTER  •  (*)  TRANS 


Dimensions 


DIMENSION 

DIMENSION 


XP  (3.3) 
XN  (3,3) 


DATA  STATEMENTS 


.ROTATE  (3,3) 


DATA  DEGREE 


DATA  SINEPS 


Degrees  to  radians. 


/  1.7453  29251  99432  96  D  -2  / 


Sine  of  obliquity  of  ecliptic. 


/  3.978  3951  D  -1 


Constant  part  of  nutation  matrix. 


DATA  XN  (  1  ,  1  )  /  1  DO  / 

DATA  XN  (2,2)  /  1 .  DO  / 

DATA  XN  (3,3 )  /  1 .  DO  / 

DATA  XN  (  1  ,2  )  /  0  .DO  / 

DATA  XN  (2 , 1 )  /  0  .DO  / 


BEGIN  PROGRAM 


Wi'TiIW 


PRECESSION  MATRIX 


Is  precession  matrix  needed? 

IF  (TRANS  .NE.  'NUTATE')  THEN 

Compute  time  parameters. 

T  «  TIME  /  36524.2200 

TT  *  T  *  T 

Compute  precession  matrix  --  1950.0  to  requested 
Date , 


XP 

(1.1) 

* 

1  .  DO  -  ( 1 . 3D-7*T  +  2 . 9696D-4 )  *  TT 

XP 

(1,2) 

= 

((2.210-6’T  -  6.760-6)  *  T  -  2.234941D 

XP 

(1.3) 

X 

( ( 9 . 6D-7*T  ♦  2.070-6)  *  T  -  9.7169D-3) 

XP 

(2.1) 

s 

-  XP( 1 .2) 

XP 

(2.2) 

X 

1 . DO  -  ( 1 . 50 - 7 * T  +  2  4S75D-4)  *  TT 

XP 

(2.3) 

s 

-  1 . 0858D -4  *TT 

XP 

(3.1) 

X 

-  XP( 1,3) 

XP 

(3.2) 

X 

XP(  2.3) 

XP 

(3,3) 

r 

1.00  -  4.721D-5*TT 

END  IF 


NUTATION  MATRIX 


Is  nutation  matrix  needed? 

IF  (TRANS  NE,  'PRECES')  THEN 

Compute  time  parameters. 

DC  «  (TIME  -  .07700)  ♦  18262.500 
TC  *  DC  /  36525.00 

Compute  auxiliary  angles. 

0  ■  ( ( ( 4 ,6D-20*DC  *  1.5570-12)  •  DC  -  5.29539220-2)  *  DC 

*  259. 1832700)  *  DEGREE 

P  *  ( ( ( 7 . 80 - 20*0C  *  1.7000-12)  *  DC  *  2 . 6352793D+ 1  )  *  DC 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


♦  540.86823D0)  *  DEGREE 
0  =  ( (4 . 52D- 13  +  DC  +  1.971294700)  •  DC 

-  160 . 6O664D0  )  *  DEGREE 
TWOO  =  0+0 

Compute  nutation  matrix  at  requested  date. 

XN  (1.3)  *  (SIN(O)  *  (17.2327D0  +  0 . 0 1 737D0  +  TC  ) 

SIN(TWOO)  *  (0.208SD0  +  0.20-4'TC)  + 

SIN(Q)  •  (1.2729D0  +  0.13D-3+TC)  + 

SIMP  )  *  (0.2037D0  +  0 . 2D-4  +  TC ) )  •  SINEPS  /  206264. 8D0 

XN  (3,1)  •  -  XN( 1,3) 

XN  (3,2)  »  (COS(O)  *  (9.21D0  +  0.91D-3>TC) 

COS(TWOO)  •  (0.0904D0  -  0.4D-4+TC)  + 

COS(O)  *  (0.5522D0  -  0.29D-3+TC)  + 

COS(P)  •  (0.088400  -  0 . 5D-4*TC  ) )  /  206264 . 8D0 
XN  (2,3)  *  -  XN( 3,2) 

END  IF 


TRANSFORMATION  MATRIX 


Nutation. 

IF  (TRANS  . EO .  'NUTATE')  THEN 
CALL  MSET33  (ROTATE,  XN) 

Precess i on . 

ELSE  IF  (TRANS  .EO.  'PRECES')  THEN 
CALL  MSET33  (ROTATE,  XP) 

Precession  and  nutation. 

ELSE 

CALL  MAMA33  (ROTATE,  XN.XP) 

END  IF 


****** 

RETURN 

****** 


RETURN 


r 
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